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Abstract 

The  use  of  adaptive  optics  entails  the  design  of  a  controller.  This  requires  the 
development  of  a  model  of  the  plant  to  be  controlled,  which,  in  this  case,  consists  of  the 
atmosphere  through  which  light  is  traveling.  In  optics,  Zernike  polynomials  are  used  as  a 
basis  set  for  the  expansion  of  wavefront  phase  distortions.  Due  to  the  turbulence  induced 
stochastic  nature  of  the  underlying  process  involved,  the  spatial-temporal  correlation 
functions  of  the  Zernike  polynomial  phase  expansion  coefficients  must  be  evaluated  if  a 
proper  stochastic  model  of  the  plant  is  to  be  developed  and  adaptive  optics  is  to  be 
employed.  In  this  work,  these  correlation  functions  are  developed  using  a  layered 
atmospheric  model  which  takes  into  account  wind  effects  and  anisoplanatism. 
Calculations  are  provided  for  the  first  few  low-order  Zernike  modes.  Using  these 
correlation  functions,  an  underlying  linear,  stochastic,  dynamical  system,  which 
represents  the  atmosphere  and  is  adequate  for  control  synthesis,  is  identified.  Within  an 
acceptable  error  bound,  the  correlation  functions  of  this  system  are  representative  of  the 
calculated  functions.  The  deformable  mirror  is  also  modeled,  output  equations  are 
specified,  and  the  complete  system  is  constructed.  This  system,  in  turn,  provides  the 
basis  for  the  employment  of  advanced  control  and  estimation  concepts.  The  control 
objective  is  to  apply  the  estimated  conjugate  phase  to  the  deformable  mirror  so  that,  at  the 
target,  the  outbound  wavefront  distortion  is  minimized  and  the  Strehl  ratio  is  maximized. 
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ALGORITHM  DEVELOPMENT  FOR  ON-LINE 
CONTROL  OF  THE  AIRBORNE  LASER 


/.  Introduction 


1.1  Overview 

For  years,  scientists  have  been  concerned  about  the  effects  of  atmospheric 
turbulence  on  the  propagation  of  light.  Atmospheric  turbulence  has  the  effect  of 
randomly  varying  the  index  of  refraction  of  the  propagation  medium.  These  fluctuations 
cause  the  familiar  “twinkling”  of  stars  and  “shimmering  of  heat”  over  hot  pavement  [25]. 
Due  to  the  randomly  varying  index  of  refraction,  atmospheric  turbulence  causes  random 
wavefront  distortion  which  in  turn  degrades  the  performance  of  optical  imaging  systems, 
below  the  limit  set  by  diffraction  theory.  However,  it  has  been  shown  that  real-time 
wavefront  compensation,  commonly  called  adaptive  optics,  can  improve  system 
performance  [11].  Adaptive  Optics  (AO)  is  a  broad  term  used  to  cover  a  system  which, 
in  real  time,  measures  the  optical  aberrations  induced  in  the  wavefront  over  the 
propagation  path  and  corrects  for  their  adverse  effects. 

The  root  of  all  evil  in  the  propagation  of  light  in  the  atmosphere  is  the  existence  of 
fluctuations  in  the  index  of  refraction.  Index  of  refraction  fluctuations  are  due  to 
turbulent  air  motion  which  obtains  its  energy  from  differential  heating  and  cooling  of  the 
atmosphere.  These  fluctuations  vary  both  spatially  and  temporally.  Because  of  this, 
when  an  optical  wave  propagates  through  a  random  medium  (like  the  atmosphere),  the 
wave  undergoes  chaotic  changes  in  its  phase,  amplitude,  and  other  parameters  [28]. 


Stochastic  models  for  turbulence  specify  second-order  statistics,  such  as  the  spatial 
correlation,  of  the  index  of  refraction  fluctuations.  Hence,  it  follows  that  the  turbulence- 
induced  wavefront  phase  and  amplitude  deformations  can  be  characterized  by  their  spatial 
correlation  properties.  Therefore,  the  statistical  properties  of  the  phase  and  amplitude 
fluctuations  must  be  determined. 

Adaptive  optics  is  a  method  of  compensating  for  the  random  phase  and  amplitude 
fluctuations  of  an  optical  wave  propagating  through  turbulence.  An  adaptive  optics 
system  consists  of  three  basic  components: 

1.  Wavefront  Sensor  (WFS)  which  is  a  device  to  measure  wavefront 
deformations. 

2.  Wavefront  modifying  device  which  corrects  wavefront  deformations, 
commonly  called  a  Deformable  Mirror  (DM). 

3.  Controller  which  processes  measurements  and  provides  commands  to  the 
wavefront  modifying  device. 

In  general  terms,  adaptive  optics  involves  the  measurement  and  control  of  wavefronts  in 
real  time  in  order  to  deposit  light  effectively  on  a  detector  or  target. 

In  the  most  basic  current  application,  adaptive  optics  is  used  in  ground-based 
imaging  systems.  In  this  case,  the  receiving  telescope  is  fixed  on  the  ground  and  is 
directed  to  an  astronomical  object.  This  scenario  is  depicted  in  Figure  1. 
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Figure  1 :  Conventional  adaptive  optics  system  used  in  astronomy. 


In  Figure  1 ,  the  three  main  components  of  an  adaptive  optics  system  are  shown. 
Specifically,  these  are  a  wavefront  modifying  device  or,  in  this  case,  a  deformable  mirror 
(the  actuator),  a  wavefront  sensor,  and  an  information  processing  device  (controller).  The 
WFS  accepts  light  and  provides  a  measurement  of  the  wavefront  slope,  while  the 
controller  accepts  the  output  of  the  WFS  and  converts  it  to  appropriate  control  signals 
sent  to  the  deformable  mirror.  Indeed,  an  AO  enhanced  telescope  is  a  feedback  control 
system  in  which  the  (uncertain)  plant  is  the  turbulent  atmosphere.  Atmospheric 
turbulence  induces  phase  shifts  in  the  light’s  wavefront  which  are  measured  by  the  WFS, 
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processed  by  the  controller,  and  corrected  by  the  DM.  Hence,  the  feedback  control 
system  is  a  regulator  -  the  control  objective  being  to  drive  the  atmospheric  turbulence- 
induced  phase  distortion  to  zero.  In  the  scenario  depicted  in  Figure  1 ,  the  WFS  accepts 
light  reflected  from  the  deformable  mirror.  In  this  case,  the  WFS  measures  the  wavefront 
distortion  in  the  corrected  wavefront.  Thus,  the  WFS  directly  provides  a  measurement  of 
the  error  signal,  as  is  required  in  conventional  SISO  control  systems. 

Figure  2  shows  a  modified  adaptive  optics  system  that  can  be  used  in  astronomy. 


TURBULENCE 


Figure  2:  Modified  adaptive  optics  system  used  in  astronomy. 


In  this  case,  the  WFS  measures  the  distortion  in  the  incident  light’s  wavefront.  The 
incident  radiation  may  be  either  the  object  being  imaged  or  a  beacon.  The  beacon  is  an 
object,  either  a  natural  star  or  a  laser  guide  star,  which  provides  light  to  the  WFS.  Many 
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times  beacons  are  necessary  because  the  object  being  imaged  does  not  provide  an 
adequate  amount  of  light.  The  metrological  laser  in  Figure  2  is  used  to  provide  a  minor 
feedback  loop  which  is  useful  when  the  dynamics  of  the  deformable  mirror  are  not  known 
precisely.  Figures  3  and  4  show  the  AO  control  systems  for  the  configurations  shown  in 
Figures  1  and  2. 


Figure  3.  Control  system  which  corresponds  to  the  AO  arrangement  shown  in  Figure  1. 


Figure  4.  Control  system  which  corresponds  to  the  AO  arrangement  shown  in  Figure  2. 
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In  Figures  3  and  4,  the  controller  is  denoted  by  G;  the  atmospheric  model  is  denoted  by 
PLANT;  w  denotes  white  process  noise  which  drives  the  plant;  the  phase  distortion  of  the 
light  deposited  on  the  detector  is  y;  and  v,  vj,  and  V2  are  the  WFS  white  noise  signals. 

The  subscript  ‘c’  denotes  the  command  signal,  which  is  zero  in  the  regulator  control 
systems  used  in  AO. 

If  the  AO  components  discussed  above  worked  perfectly,  the  AO  system  would 
achieve  ideal  performance,  namely,  the  diffraction  limit.  However,  there  are  numerous 
factors  which  degrade  the  performance  of  an  AO  system.  The  fundamental  limiting 
factors  include  [23]: 

1.  Finite  light  levels  in  the  WFS. 

2.  Finite  degrees  of  freedom  of  the  DM. 

3.  Finite  spatial  sampling  of  the  WFS. 

4.  Finite  temporal  response  of  the  AO  system. 

5.  Anisoplanatism. 

Until  recently,  factors  (1)  -  (4)  were  severe  problems.  Fortunately,  with  advances  in 
WFS,  DM,  and  AO  hardware  design,  these  degradations  play  a  much  smaller  role.  The 
same  cannot  be  said  for  anisoplanatism.  Anisoplanatism  is  a  term  referring  to  the  object 
and  beacon  wavefronts  sampling  different  volumes  of  the  atmosphere.  Anisoplanatism 
results  when  the  object  and  beacon  have  a  spatial  separation.  This  spatial  separation 
between  the  beacon,  or  reference,  and  object  causes  the  optical  paths  from  the  reference 
and  the  object  to  traverse  different  regions  of  the  atmosphere,  resulting  in  distinct 
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wavefront  perturbations  for  each  wavefront  [23].  Figure  5  displays  a  common 
anisoplanatic  scenario. 


BEACON  OBJETCT 


Figure  5:  Anisoplanatism. 


However,  just  because  a  spatial  separation  exists  between  the  reference  and  object 
does  not  mean  that  anisoplanatic  effects  are  observed.  In  fact,  the  two  objects  may  lie 
within  the  isoplanatic  angle,  defined  as  the  maximum  angular  separation  for  which  the 
turbulence-induced  wavefront  deformations  for  the  object  and  beacon  wavefronts  are  still 
essentially  the  same.  In  the  visible  spectrum,  the  object  being  imaged  must  be  within  5  to 
10  (irad  of  the  beacon  [23],  At  infrared  wavelengths,  the  isoplanatic  angle  is  significantly 
larger,  on  the  order  of  100’s  of  |0.rad. 
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1.2  AirBorne  Laser 


The  AirBorne  Laser  (ABL)  is  an  atmospherically  based  directed  energy  weapon 
system  currently  being  developed  by  the  United  States  Air  Force.  This  program  uses  a 
high-energy  laser,  mounted  in  an  airplane,  to  deposit  energy  on  and  disable  a  theater 
ballistic  missile  in  boost  phase.  The  ABL  scenario  is  significantly  different  from  typical 
optical  imaging  systems.  First,  the  ABL  operates  at  an  altitude  of  about  40,000  ft.  Unlike 
ground-based  imaging  systems,  in  which  the  wavefronts  propagate  perpendicularly  to  the 
earth’s  surface,  the  propagation  paths  of  the  beacon  and  High-Energy  Laser  (HEL)  are 
approximately  horizontal  to  the  earth’s  surface.  Second,  there  is  an  enormous  amount  of 
motion  in  this  problem.  Both  the  aperture  (mounted  in  the  airplane)  and  target  are 
moving  at  large  velocities.  Hence,  unlike  conventional  imaging  systems  in  which  the 
wind  introduces  a  temporal  dependence  in  the  index  of  refraction  fluctuations,  dynamic 
aperture  and  target  motion  account  for  the  majority  of  temporal  effects. 

The  arrangements  shown  in  Figures  1  and  2  do  not  allow  for  anisoplanatism  and 
do  not  account  for  time  delays  inherent  in  the  control  system.  For  use  in  the  ABL,  the 
adaptive  optics  configurations  shown  in  Figures  1  and  2  must  be  modified.  Figure  6 
shows  the  ABL’s  adaptive  optics  arrangement.  In  this  situation,  there  are  two  wavefronts 
of  interest.  First  is  the  beacon,  incident,  or  inbound,  radiation.  The  beacon  wavefront  is  a 
spherical  wavefront  emanating  from  the  sharp  tip  of  the  missile  -  or  a  fixed  anchor  point. 
This  is  unlike  the  astronomical  case  in  which  radiation  from  the  laser  guide  star  beacon  is 
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Figure  6:  ABL  adaptive  optics  system. 


scattered  from  the  layer  in  the  ionosphere.  Hence,  in  astronomy,  the  beacon  is  subject  to 
constant  turbulence-induced  fluctuations  due  to  the  upward  passage  of  the  beacon’s  laser 
beam  through  the  turbulent  atmosphere.  This  does  not  allow  for  tilt  estimation  (and 
correction)  in  the  astronomical  application  of  straight-forward  AO.  This  is  not  so  in  the 
ABL  case  in  which  the  target’s  tip  is  flooded  by  laser  illumination  and  most  of  the 
radiation  is  reflected  off  the  fixed  tip  of  the  missile.  The  beacon  wavefront  is  the  same  as 
the  ones  displayed  in  Figures  1  and  2.  Second  is  the  high-energy  beam  which  is  used  to 
deposit  the  corrected,  or  controlled,  high-energy  laser  on  the  target.  As  in  Figure  2,  the 
metrological  laser  is  used  to  provide  a  minor  feedback  loop  if  the  mirror  dynamics  are  not 
known  precisely. 

Atmospheric  turbulence  will  induce  random  temporal  variations  in  the  laser 
beam’s  wavefront  which,  in  turn,  will  reduce  the  light  intensity  at  the  target  and  hence  the 


9 


directed  energy  weapon’s  effectiveness.  This  fact  provides  the  motivation  for  the 
employment  of  AO  in  order  to  correct  for  the  deleterious  effects  caused  by  atmospheric 
turbulence. 

In  order  to  apply  AO  to  the  AirBorne  Laser  (ABL),  a  model  of  the  plant  to  be 
controlled  must  be  developed.  The  plant,  in  this  case,  is  the  turbulent  atmosphere  through 
which  light  is  traveling.  Because  of  atmospheric  turbulence,  the  plant,  which  models  the 
atmosphere,  will  be  stochastic. 

1.3  Historical  Overview 

Pioneering  work  in  optical  wave  propagation  through  atmospheric  turbulence  was 
performed  in  the  1950’s  by  Tatarski  [27].  In  his  work,  Tatarski  considered  the  medium 
between  the  source  and  aperture  as  a  continuum.  However,  the  most  common  forms  of 
the  Power  Spectral  Density  (PSD)  of  the  index  of  refraction  fluctuations  contain 
singularities  at  the  origin.  Thus,  evaluation  of  the  necessary  integrals  requires  special 
care.  In  the  simple  astronomical  scenario  considered  by  Tatarski,  integrations  in  the  z- 
direction,  i.e.,  the  direction  of  wavefront  propagation,  could  be  performed  in  closed-form. 
In  the  process,  the  singularity  in  the  PSD  of  the  index  of  refraction  fluctuations  did  not 
pose  any  difficulty;  however,  no  such  analytical  techniques  present  themselves  in  the 
more  dynamic  ABL  scenario.  A  layered  (discrete)  atmospheric  model  can  be  used  to 
avoid  integration  in  the  direction  of  wavefront  propagation.  In  this  model,  the 
atmosphere  is  assumed  to  be  composed  of  a  finite  number  of  properly  separated  layers  so 
that  the  spatial  correlation  can  be  described  as  the  sum  of  the  contributions  from  all 
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layers.  Hence,  the  layered  atmospheric  model  requires  simpler,  two-dimensional 
integrations  as  opposed  to  the  three-dimensional  integrations  in  Tatarski’s  original  work. 
Other  seminal  work  in  optical  propagation  through  atmospheric  turbulence  is  due  to 
Kolmogorov  [14]  and  Fried  [4,  5].  Kolmogorov  developed  a  statistical  model  for  the 
structure  of  turbulent  air  flows.  His  model  is  widely  used  as  the  spatial  PSD  of  the  index 
of  refraction  fluctuations.  Fried  has  extended  the  results  of  Tatarski  and  also  developed 
the  seeing  cell  size,  commonly  called  the  Fried  parameter,  r0.  The  seeing  cell  size  is 
interpreted  as  the  aperture  size  beyond  which  further  increases  in  diameter  result  in  no 
further  increases  in  resolution. 

Zernike  polynomials  have  been  extensively  used  for  analysis  in  optical  systems  [6, 
12,  18].  These  orthonormal  polynomials,  which  are  defined  within  a  circle  of  unit  radius, 
form  a  complete  set  and  represent  many  specific  phase  aberrations  [23].  Due  to  the  fact 
that  the  ABL  application  requires  the  consideration  of  two  distinct  wavefronts  at  two 
different  time  instants,  the  spatial-temporal  correlation  of  the  Zernike  polynomials’  phase 
expansion  coefficients  of  the  wavefronts  of  interest  is  required.  For  implementation  in 
adaptive  optics  systems,  some  partial  results  for  these  correlation  functions  have  been 
provided  by  Takato  and  Yamaguchi  [26],  Valley  and  Wandzura  [30],  Whiteley,  et.  al 
[35],  and  Whiteley  [36], 
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1.4  Problem  Statement 


The  application  of  AO  to  the  ABL  is  a  vast  problem.  In  no  way  could  one 

dissertation  attempt  to  tackle  the  entire  problem.  Therefore,  a  viable,  but  challenging 

problem  statement  for  this  dissertation  is  as  follows: 

The  objectives  are  to  develop  fully  a  unified  approach  to  beam  control  for  the 
ABL  so  that  advanced  control  techniques  can  be  applied,  and  to  quantify  the 
performance  improvements  that  can  be  achieved  by  utilizing  these  techniques. 

To  this  end,  a  generic  ABL  engagement  geometry  is  selected  as  the  basis  for  this  problem. 

Simplifications  can  be  made  to  transform  this  work  into  special  cases  which  are 

applicable  to  imaging  and  laser  communications.  For  this  general  geometry,  the 

performance  degradation  due  to  turbulence-induced  optical  aberrations  is  investigated  and 

the  associated  optimal  estimation  and  control  theory  is  developed.  The  results  obtained  in 

this  research  are  for  a  layered,  homogeneous,  isotropic,  wide-sense  stationary  turbulence 

model.  For  a  homogeneous,  isotropic  random  process,  the  autocorrelation  function  of 

X(r ) ,  denoted  Tx  (?) ,  is  a  function  only  of  the  scalar  distance  r  =  |r| .  Tilt  compensation 

is  the  main  concern  in  this  problem,  with  model  inadequacies  and  time  delays  also 
included.  The  control  objective  is  to  provide  commands  to  the  DM  so  that  the  estimated 
conjugate  phase  can  be  applied  to  the  mirror  such  that,  at  the  target,  the  high-energy 
laser’s  wavefront  distortion  is  minimized  and  the  Strehl  ratio  is  maximized.  In  this  way, 
the  HEL  is  correctly  pointed  to  the  aim  point.  The  Strehl  ratio  is  a  performance  metric 
defined  as  the  irradiance  with  wavefront  aberrations  present  divided  by  the  diffraction- 
limited  irradiance  on  the  optical  axis.  A  perfect  Strehl  ratio  is  one,  while  typical  imaging 
systems  have  Strehl  ratios  well  below  the  theoretical  limit;  however,  an  exact  number  is 
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difficult  to  present  because  many  factors  contribute  to  the  Strehl  ratio.  A  good  Strehl 

ratio  would  be  0.6  or  higher. 

1.5  Key  Results 

The  analysis  conducted  in  this  research  has  led  to  the  following  key  results: 

•  A  unified  framework  for  the  application  of  adaptive  optics  to  the  ABL. 

•  A  fully  developed  engagement  geometry  for  the  ABL  encompassing  target  and 
aperture  motion  as  well  as  wind  effects. 

•  A  modified  Frozen  Flow  Hypothesis  (FFH)  which  takes  into  account  translation  and 
rotation  of  the  aperture  and  target. 

•  A  new  approach  to  developing  an  atmospheric  model  which  directly  uses  correlation 
kernel  data  instead  of  designing  shaping  filters  for  the  atmospheric  states.  This  results 
in  a  lower-order,  and  therefore,  more  practical  atmospheric  model. 

•  Development  of  an  atmospheric  model  that  provides  a  good  representation  of  the 
atmosphere.  This  is  validated  by  examining  the  robustness  of  the  Kalman  filter  based 
upon  that  model. 

•  A  complete  AO  control  system  which  accounts  for  anisoplanatism  and  time  delays. 

•  Significant  improvement  in  the  reduction  of  wavefront  phase  deformations  that  can  be 
obtained  by  simply  performing  tilt  correction.  This  result  is  displayed  using 
simulations  of  the  entire  ABL  AO  control  system. 
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•  Strehl  ratios  that  are  improved  by  about  0.15  using  tilt  compensation  as  compared  to 
the  no  compensation  case.  In  terms  of  simple  conjugation,  Strehl  ratios  are  improved 
by  about  0.04. 

1.6  Organization 

Due  to  the  atmospheric  turbulence-induced  stochastic  nature  of  the  underlying 
processes  involved,  the  spatial-temporal  correlation  functions  of  the  Zernike  polynomial 
phase  expansion  coefficients  must  be  evaluated  if  a  proper  stochastic  model  of  the  plant  is 
to  be  developed  and  adaptive  optics  is  to  be  employed.  In  this  research,  these  correlation 
functions  are  developed  using  a  layered  atmospheric  model,  and  calculations  for  the  first 
few  low-order  Zernike  modes  are  performed.  The  spatial-temporal  correlation  of  the 
Zernike  phase  expansion  coefficients  is  developed  for  a  generic  adaptive  optics 
configuration.  In  fact,  the  general  situation,  as  it  applies  to  the  ABL  scenario,  is 
considered.  Specifically,  two  separate  wavefronts  at  two  distinct  time  instants  are 
considered.  With  the  dynamic  scenario  in  mind,  a  novel  modified  frozen  flow  hypothesis 
which  accounts  for  translation  as  well  as  rotation  of  the  aperture  and  target  is  also 
developed. 

The  use  of  AO  entails  the  design  of  a  controller.  This  requires  the  development  of 
a  model  of  the  plant  to  be  controlled.  In  AO,  the  plant  consists  of  the  atmosphere  through 
which  light  is  traveling.  Moreover,  a  distinct  feature  of  the  AO  control  application  is  the 
presence  of  random  signals  in  the  plant.  Thus,  in  order  to  apply  control  concepts  to  the 
ABL  scenario  and  realize  the  benefits  of  AO,  an  underlying  linear,  stochastic,  dynamical 
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system  model  which  generates  the  correlation  functions  of  the  Zernike  expansion 
coefficients  and  is  adequate  for  control  synthesis  must  be  developed.  This  work 
addresses  the  attendant  realization  problem.  Using  the  calculated  correlation  functions, 
an  atmospheric  model  (plant)  is  constructed.  This  model  will  then  serve  in  a  Kalman 
filter-based  optimal  stochastic  controller  (Linear  Quadratic  Gaussian  controller) 
mechanization  of  an  adaptive  optics  control  system  for  the  ABL.  In  the  ABL  paradigm, 
the  deformable  mirror  is  the  actuator  element  for  the  control  system.  Thus,  a  model  of 
the  DM  compatible  with  the  modal  approach  is  developed.  An  equivalent  discrete-time 
model  of  the  realized  continuous-time  system  is  constructed  and  the  complete  adaptive 
optics  control  system  is  discussed  and  simulated. 

This  dissertation  is  organized  as  follows.  Atmospheric  models  and  Zernike 
polynomials  are  discussed  in  Chapter  2.  In  Chapter  3,  the  ABL-specific  engagement 
scenario  is  introduced,  including  all  of  the  necessary  geometrical  relations.  Next,  the 
attendant  auto-  and  cross-correlations  of  the  atmosphere’s  Zernike  polynomial  phase 
expansion  coefficients  are  derived  in  Chapter  4.  Also,  the  correlation  functions  are 
calculated  and  a  linear,  stochastic,  dynamical  system  is  synthesized  using  the  developed 
correlation  functions.  Chapter  5  contains  a  description  of  the  modeling  of  the  deformable 
mirror  dynamics,  using  a  modal  approach.  The  complete  system  model  as  well  as  the 
equivalent  discrete-time  system  are  also  discussed  in  Chapter  5.  Chapter  6  describes  the 
Kalman  filter  and  Linear  Quadratic  Gaussian  (LQG)  controller  design  while  Chapter  7 
provides  results  of  the  control  system  simulation  as  well  as  a  discussion  of  these  results. 
Last,  appendices,  a  succinct  list  of  relevant  references,  and  a  vita  are  provided. 
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1.7  Summary 

In  this  chapter,  the  compensation  of  atmospheric  turbulence  using  adaptive  optics 
has  been  introduced.  The  main  components  of  an  AO  system  as  well  as  some  common 
AO  configurations  have  also  been  considered.  Some  of  the  most  prevalent  performance 
degrading  factors,  with  emphasis  on  anisoplanatism,  have  also  been  discussed.  To 
achieve  the  goal  of  successfully  applying  AO  to  the  ABL,  a  general  engagement  geometry 
must  be  considered.  This  allows  development  of  the  correlation  between  the  inbound 
(beacon)  and  outbound  (HEL)  wavefronts.  Advanced  control  concepts  can  then  be  used 
to  compensate  for  optical  wavefront  aberrations. 

In  the  next  chapter,  second-order  statistics  of  the  index  of  refraction  fluctuations 
will  be  discussed.  In  particular,  the  two  most  common  spectral  representations  of  the 
fluctuations,  namely,  the  Kolmogorov  and  von  Karman  spectra,  will  be  delineated.  The 
layered  atmospheric  model  and  Zernike  polynomials  will  also  be  covered. 
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II.  Atmospheric  Models  and  Zemike  Polynomials 


2.1  Introduction 

Atmospheric  models  and  Zemike  polynomials  will  be  discussed  in  the  upcoming 
sections.  In  order  to  describe  the  atmosphere  mathematically,  the  physical  nature  of 
turbulence  must  be  understood.  In  particular,  index  of  refraction  fluctuations,  the 
Kolmogorov  spectaim,  and  the  von  Karman  spectrum  will  be  covered.  Also,  the  layered 
atmospheric  model  will  be  examined.  Following  this,  a  description  of  Zernike 
polynomials  and  their  application  to  optics  will  be  included. 


2.2  Index  of  Refraction  Fluctuations 

Index  of  refraction  fluctuations  arise  from  turbulent  air  motion  which  obtains  its 
energy  from  differential  heating  and  cooling  of  the  atmosphere.  Turbulent  eddies  are 
pockets  of  air  which  have  a  uniform  index  of  refraction.  The  statistical  description  of  the 
number  and  size  of  these  eddies  is  given  by  the  Power  Spectral  Density  (PSD)  of  the 
randomly  fluctuating  index  of  refraction  term,  n,  (r,t),  where  r  is  a  three-dimensional 
position  vector  and  t  is  time.  Denoting  this  PSD  by  Wn|  (k,Q  (explicitly  a  function  of 


spatial  wavenumber  K  and  altitude  £,),  two  widely  used  forms  can  be  delineated,  namely, 
the  Kolmogorov  and  von  Karman  spectra  [13]: 

/  s') 


V3  r 


wnK(K,CT 


2(2n 


C„(Qk  3 


(1) 


17 


and 


where  Kq  =  27t  /  L0,  L0  is  the  outer  scale  of  turbulence,  Cn2(Q  is  structure  constant  of  the 
index  of  refraction  fluctuations,  r(«)  is  the  gamma  function,  and  the  superscripts  ‘K’  and 

‘V’  represent  the  Kolmogorov  and  von  Karman  spectra,  respectively.  Cn2(Q  can  be 
thought  of  as  a  term  which  characterizes  the  strength  of  the  index  of  refraction 
fluctuations.  Many  Cn2(Q  profiles  exist,  for  example,  the  SLC-Day,  Hufnagel- Valley, 
and  Greenwood  models  [10,  21].  The  L0  term  represents  the  characteristic  dimension  of 
the  largest  turbulent  eddies  which  break  up  following  Kolmogorov  theory.  The  von 
Karman  spectrum  will  be  used  extensively  in  the  upcoming  work  because  it  does  not 
contain  the  non-integrable  pole  at  k  =  0  which  is  present  in  the  Kolmogorov  spectrum. 

Thus  far,  only  the  spectral  characteristics  of  the  index  of  refraction  fluctuations 
have  been  considered.  These  spectral  representations  have  corresponding  spatial 
representations,  with  the  two  being  related  through  the  Fourier  transform.  As  already 
stated,  the  second-order  statistics  of  the  Zemike  polynomial  phase  expansion  coefficients 
are  of  interest.  However,  to  obtain  this  expression  requires  evaluation  of  the  second- 
order  statistics  of  the  index  of  refraction  fluctuations.  In  order  to  examine  this,  two 
definitions  related  to  three-dimensional  random  processes  will  be  required.  First,  let  K(r) 
be  a  three-dimensional  random  process  where  r  is  a  position  vector.  The  process  K(  r)  is 
said  to  be  homogeneous  if  the  autocorrelation  of  K(r), 
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rK(f„r2)  =  E{K(f,)K(?2)},  (3) 

where  E{  • }  denotes  the  expectation  operator,  is  a  function  only  of  the  difference  r,  -  r2. 
Second,  a  random  process  is  isotropic  if  its  autocorrelation  function  is  spherically 
symmetric.  In  other  words,  for  a  homogeneous,  isotropic  random  process,  the 
autocorrelation  function,  rK(r) , is  a  function  only  of  the  scalar  distance  r  =  |r|,  where  |«| 

denotes  the  magnitude  of  the  corresponding  vector.  In  all  cases,  the  index  of  refraction 
fluctuations  in  the  atmosphere  are  assumed  to  be  homogeneous  and  isotropic. 

2.3  Layered  Atmospheric  Model 

In  this  model,  the  atmosphere  is  assumed  to  be  composed  of  a  finite  number  of 
properly  separated  layers.  The  term  properly  refers  to  the  layers  being  treated  as 
independent.  A  more  thorough  discussion  of  independent  layers  is  presented  later  in  this 
section.  Basically,  the  layered  model  can  be  thought  of  as  the  discrete  counterpart  to  the 
continuous  model.  Hence,  integrations  along  the  propagation  path  (z-direction)  in  the 
continuous  model  are  replaced  by  summations  in  the  discrete  model.  Therefore,  the 
desired  spatial  correlation  between  wavefronts  can  be  described  as  the  sum  of  the 
contributions  from  all  layers.  This  layering  of  the  turbulence  significantly  simplifies  the 
calculations  to  follow. 

For  a  layered  model,  the  Cn2(C)  profile  is  broken  up  into  a  finite  number  of  slabs, 
with  each  slab  characterized  by  a  turbulence  strength  which  is  nearly  constant  in  each 
layer.  Letting  Cn2(Ci)>  Ci>  an<3  A£j  be  the  structure  constant,  altitude,  and  thickness  of  the 
ith  turbulent  layer,  the  structure  constant  and  altitude  are  selected  in  such  a  way  that  the 
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zeroth  through  seventh  moments,  of  the  continuous  model  match  the  discrete  model 
[23]: 


J?'c;©<ic=xt?c;(c,)AC,  •  (« 

0  'l 

In  eq.  (4),  L  is  the  propagation  path  length  and  0  <  t,  <  7.  Moments  zero  through  seven 
are  used  to  provide  enough  fidelity  in  the  layered  model.  Higher  order  moments  could  be 
used,  but  the  rate  of  return  is  negligible.  The  weights  and  heights  of  the  most  common 
four-layer  models  have  been  calculated  by  Troxel,  et.  al.  [29].  Also,  it  should  be 
mentioned  that  A£j  surrounds  C,,  symmetrically. 

If  diffraction  effects  are  assumed  to  be  negligible,  then  the  turbulence-induced 
phase  at  position  x  can  be  written  as 


u  ' 

where  X  is  the  wavelength.  Equation  (5)  displays  the  computation  of  wavefront  phase 
when  the  atmosphere  is  considered  a  continuum.  Since  <b(x)  is  the  sum  of  a  large 
number  of  zero-mean  random  variables,  namely  n,  (x,  z) ,  then,  by  the  Central  Limit 


(5) 


Theorem  [3],  <E>(x)  will  be  zero-mean  and  Gaussian. 

For  a  layered  atmospheric  model,  the  expression  analogous  to  eq.  (5)  is 

4>(x)=^5>Cin,(x,l;i).  (6) 

A  i 

By  the  Fourier  transform  relation  between  the  PSD  and  autocorrelation,  the  PSD  of  the 
turbulence-induced  phase  is 
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W0(k)  =  FT[e{<F(x)<I>(x  +  Ax)}] 


=  S  A^i  (  y)?  a^'  (y)ft[e^  n  1  frCi  )n>  (* + M’  )}] 


SA<f)w„(K,,) 


(7) 


where  FT{*}  denotes  the  Fourier  transform.  The  term  in  the  summation  is  the  PSD  of  the 
phase  in  the  i,h  layer,  denoted  W^k,^).  Thus, 


w4K)=Iwt(K,g. 

i 

Substituting  eq.  (2)  into  eq.  (7)  yields  the  von  Karman  spectrum  for  the  phase  in  the  i 
layer: 
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For  £  =  0,  eq.  (4)  implies  that 
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The  term  in  eq.  (10)  can  be  identified  as  a  turbulence-dependent  layer  weight  Wj: 
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Defining  the  atmospheric  coherence  diameter,  or  Fried  parameter  [4],  r0,  as 
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and  using  this  result  and  eq.  (1 1)  in  eq.  (9)  gives 
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Equation  (13)  provides  the  final  form  of  the  von  Karman  spectrum  and  will  be  used  in 
subsequent  calculations. 

When  using  a  layered  atmospheric  model,  the  problem  of  calculating  the 
correlation  properties  of  the  field  perturbations  reduces  to  calculating  those  for  a  single 
layer  and  then  extending  the  results  to  account  for  all  layers.  Tatarski  [27],  Goodman  [8], 
and  Troxel,  et.  al.  [29]  have  argued  that  the  layers  can  be  treated  as  approximately 
independent  if  the  separation  of  layer  centers  is  greater  than  the  largest  distance  between 
field  points  in  the  pupil  [23].  This  allows  the  spatial  correlation  of  phase  perturbations  to 
be  calculated  for  each  layer  separately,  with  the  results  being  added  to  account  for 
propagation  through  the  entire  turbulent  region.  The  independence  condition  will  be 
satisfied  by  choosing  the  locations  of  the  turbulent  layers  such  that  the  separation  of  layer 
centers  is  greater  than  the  largest  distance  between  field  points  in  the  pupil,  i.e.,  the 
aperture  diameter  in  this  case. 

The  results  in  this  dissertation  are  derived  using  geometrical  optics  which  is 
equivalent  to  assuming  that  near  field  conditions  exist.  Near  field  conditions  exist  if  the 
total  thickness  of  the  turbulent  region  satisfies  [37] 
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where  N  is  the  number  of  layers,  X  is  the  wavelength,  and  A  xmax  is  the  maximum 
separation  in  the  aperture.  When  near  field  conditions  exist,  the  turbulence  region  is 
close  enough  to  the  aperture  plane  that  refraction  caused  by  index  of  refraction  variations 
can  be  neglected  [23],  This  allows  the  correlation  between  the  inbound  and  outbound 
wavefronts  to  be  derived  based  on  straight  ray  path  calculations.  As  the  turbulent  region 
extends  outward  from  the  aperture  plane,  refraction  cannot  be  ignored.  In  this  case,  the 
propagation  paths  cross  and  interfere  with  each  other,  which  in  turn  causes  amplitude 
perturbations.  This  scenario  is  referred  to  as  far  field  turbulence.  In  the  case  of  near  field 
conditions,  the  perturbations  in  the  aperture  plane  are  entirely  phase  perturbations.  Near 
field  conditions  will  exist  in  the  ABL  scenario,  which  implies  that  only  phase 
perturbations  will  be  investigated. 

2.4  Zernike  Polynomials 

Zernike  polynomials,  denoted  by  Z-  (x) ,  are  a  complete  set  of  orthonormal 
polynomials  defined  within  a  circle  of  unit  radius.  They  are  commonly  used  in  optics 
because  they  represent  many  specific  phase  aberrations  [1].  In  the  study  of  optics,  it  is 
necessary  to  have  some  means  to  express  the  phase  distortion  in  a  wavefront. 

Considering  a  circular  aperture  with  an  incident  wavefront  propagating  through 
turbulence,  the  absolute  phase,  <h(x,  y,  t) ,  of  the  wave  within  the  aperture  is  a  function  of 
time  and  the  spatial  coordinates  within  the  aperture.  Here,  x  and  y  are  normalized  x- 


and  y-coordinates,  for  describing  the  plane  of  the  aperture,  such  that 
0  <  |x|  <  1 ,0  <  |y|  <  1 .  Because  Zernike  polynomials  are  orthonormal  only  on  a  unit 

circle,  the  aperture’s  radius  must  be  normalized  to  apply  Zernikes  to  non-unit  radius 
apertures.  The  appropriate  transformation  is 


(15) 


where  R  is  the  radius  of  the  aperture. 

Since  Zernike  functions  are  a  complete  set  which  span  the  functional  space 
containing  3>(x,y,t) ,  the  phase  can  be  written  as  a  linear  combination  of  Zernikes,  that 


is, 

<f)(x,y,t)  =  0>(Rx,Ry,t)  =  ^a,  (t)Zj  (x,y)  (16) 


where  aj(t)  are  the  Zernike  polynomial  phase  expansion  coefficients.  Letting  (Rx,  Ry)  = 
Rx  and  (x,  y)  =  x ,  eq.  (16)  can  be  equivalently  written  as 

3>(Rx,  t)  =  ^  aj  (t)Zi  (x) .  (17) 

i 

Other  basis  functions  exist,  for  example,  Legendre  [24]  and  Karhunen-Loeve  [18, 
34],  The  decision  to  use  Zernike  polynomials  in  this  work  is  due  simply  to  personal 
preference.  Zernike  polynomials  are  suitable  to  optics  problems  since  they  represent 
many  specific  phase  aberrations  such  as  tilt,  coma,  astigmatism,  etc.  Hence,  they  provide 
a  physical  understanding  of  the  turbulence-corrupted  wavefront.  The  drawback  of 
Zernikes  is  that  all  of  the  modes  are  not  uncorrelated  with  each  other.  Karhunen-Loeve 
basis  functions  are  simply  linear  combinations  of  Zernike  functions  and  are  uncorrelated 
[23].  However,  since  they  are  linear  combinations  of  Zernikes,  the  physical 
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understanding  of  wavefront  distortion  is  lost.  Hence,  the  decision  was  made  to  use  the 
more  common  and  physically  practical  Zemike  polynomials. 


The  expansion  coefficients  in  eqs.  (16)  and  (17)  can  be  obtained  by  projecting  the 
wavefront  phase  onto  the  space  of  Zemike  polynomials,  that  is, 

a  j  (t)  =  J  0(Rx,  t)Zj  (x)W(x)dx  (18) 

where  W(x)  is  a  weighting  function  defined  as 


if  |x|  <  1 
otherwise 


(19) 


and  the  integration  is  over  the  entire  two-dimensional  plane.  Also,  Zernike  polynomials 
form  an  orthonormal  basis  set  which  satisfies 

JW(x)Zi(x)Zj(x)dx  =  5ij  (20) 


where  8i ,  is  the  Kronecker  delta  defined  by 


8  =  j  1  * 

1J  [0  if  i*j- 


Table  1  lists  some  low-order  Zernike  polynomials,  Zj(r,0)  in  polar  coordinates,  Zj(x,y)  in 
rectangular  coordinates,  their  corresponding  radial  and  azimuthal  orders  [18],  n  and  m, 
and  the  optical  aberrations  they  represent  [1]. 
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Zemike  # 
i 

n 

m 

Zj(r,0) 

Zi(x,y) 

Aberration 

1 

0 

0 

1 

1 

Piston 

2 

■ 

1 

2r cos0 

R 

2x 

~rT 

X-tilt 

3 

■ 

1 

2r  sin0 

R 

2y 

R 

Y-tilt 

4 

2 

0 

■ 

1 

■ 

||(2x2+2y2-R2) 

Focus 

5 

2 

2 

2V6,  V 

R2  ™ 

Astigmatism 

6 

2 

2 

j§(V-yJ) 

Astigmatism 

7 

3 

1 

Vs 

J-* 

sin(0) 

^(3x2+3y2-2R2)y 

Coma 

8 

3 

■ 

Vs 

V 

-2 
/  V 

a 

cos(0) 

^-(3x2  +  3y2-2R2)x 

Coma 

Table  1.  Zernike  polynomials. 


2.5  Summary 

In  this  chapter,  index  of  refraction  fluctuations  and  the  Kolmogorov  and  von 
Karman  spectra  were  discussed.  The  layered  atmospheric  model  was  introduced  and  the 
turbulence-induced  phase,  using  this  model,  was  delineated.  Zemike  polynomials  and 
their  application  to  optics  were  also  covered. 

In  Chapter  3,  the  ABL  engagement  geometry  will  be  introduced.  The  objective  is 
to  determine  the  projections  of  the  aperture  vectors  into  the  turbulent  layers.  These 
projections  are  necessary  since  the  only  contributions  to  the  correlation  between 
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wavefronts  are  due  to  the  turbulent  slabs.  Once  the  projections  have  been  determined,  the 
separation  between  the  projected  vectors  can  be  evaluated.  The  effect  of  wind  is  also 
taken  into  account  and  the  necessary  kinematic  variables  are  delineated. 
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III.  ABL  Engagement  Geometry 


3.1  Introduction 

The  engagement  geometry  for  the  ABL  is  introduced  in  this  chapter.  Due  to  the 
dynamic  nature  of  the  ABL  engagement  scenario,  the  geometry  has  been  chosen  to  be  as 
general  as  possible.  Two  times  instants  of  interest,  t]  and  t2,  are  considered  and  both 
aperture  and  source  motion  as  well  as  wind  effects  are  included.  In  conventional  AO,  the 
latter  is  responsible  for  introducing  a  temporal  dependence  in  the  index  of  refraction 
fluctuations.  In  the  ABL  application,  dynamic  object  and  source  motion  cause  the  major 
temporal  dependence.  Since  the  object  and  source  velocities  are  much  greater  than  the 
wind  velocity,  the  wind  has  a  diminished  effect  in  this  application.  However,  wind  is  also 
included  so  that  the  general  results  derived  in  this  work  can  be  specialized  to  classical  AO 
situations  such  as  fixed  target  or  source  and  fixed  target  and  source.  Also,  a  layered 
atmospheric  turbulence  model  is  used  with  a  wind  vector  of  general  direction.  In  using  a 
layered  atmospheric  model,  it  is  necessary  to  project  the  aperture  vectors  into  the 
turbulent  layers. 

In  the  ABL,  there  are  two  wavefronts  of  interest,  namely,  the  inbound  (beacon) 
and  outbound  (HEL)  wavefronts.  Currently,  for  the  ABL,  the  inbound  wavefront’s 
wavelength  is  1.03  pm  while  the  wavelength  of  the  outbound  wavefront  is  1.05  pm. 
Therefore,  in  the  simulations  in  this  work,  it  will  be  assumed  that  the  beacon’s 
wavelength  is  1.03  pm  and  the  HEL’s  wavelength  is  1.05  pm. 
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It  is  desired  to  determine  the  correlation  between  these  two  wavefronts.  This  is 


needed  since  measurements  come  from  the  inbound  propagation  path;  that  is,  only  the 
turbulence-induced  wavefront  deformations  of  the  inbound  wavefront  are  measured.  In 
order  to  fire  the  HEL  successfully,  the  effects  of  turbulence  on  the  propagation  of  the 
outbound  wavefront  must  be  determined.  A  method  of  achieving  this  result  is  to  evaluate 
the  correlation  between  the  inbound  and  outbound  wavefronts.  Since  a  layered 
atmospheric  model  has  been  selected,  the  spatial  separation  between  inbound  and 
outbound  wavefront  points,  in  the  turbulent  layer,  is  needed.  Thus,  the  objective  in  this 
section  is  to  develop  an  expression  for  the  separation  of  projected  aperture  vectors. 
However,  to  obtain  this  quantity  it  will  be  necessary  to  project  the  aperture  vectors, 
perform  transformations  between  reference  frames,  determine  the  time  t2  kinematic 
variables  in  terms  of  their  time  ti  counterparts,  and  include  the  wind  effect.  In  this 
process,  a  novel  frozen  flow  hypothesis  needs  to  be  developed. 

3.2  ABL  Engagement 

Consider  the  simplified  ABL  engagement  geometry  depicted  in  Figure  7.  This 
figure  displays  the  aperture  and  target  at  two  time  instants  of  interest,  t|  and  t2,  where  t2  > 
ti.  Aperture  and  target  motion  are  considered  with  corresponding  constant  velocity 
vectors  VA  and  VT ,  respectively.  Constant  velocity  vectors  have  been  assumed  since  the 
sample  period,  At  =  t2  -  tj,  will  be  small  (on  the  order  of  tenths  of  milliseconds). 
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Figure  7.  Simplified  ABL  engagement  geometry. 


Two  wavefronts  are  depicted,  namely,  the  inbound  (beacon)  and  outbound  (HEL). 
Measurements  for  the  ABL  come  from  the  inbound  wavefront,  that  is,  the  tip  of  the 
missile  is  flooded  with  a  laser  and  the  returned  light  is  sampled  by  a  WFS  to  determine 
the  amount  of  distortion  in  the  inbound  wavefront’s  phase.  The  HEL  is  pointed  to  the 
aim  point,  or  that  point  where  the  beam  will  be  directed  to  disable  the  target.  It  can  be 
seen  that  there  is  a  spatial  separation  between  the  inbound  and  outbound  wavefronts, 
resulting  in  some  degree  of  anisoplanatism.  Due  to  the  fact  that  there  are  two  wavefronts, 
two  aperture  locations  at  each  time  instant  are  denoted,  namely,  RX;  and  Ru s ,  i  =  1,  2. 
Also,  it  is  assumed  that  the  aperture  accurately  tracks  the  target  at  all  times. 

The  layered  atmospheric  model  is  depicted  by  displaying  turbulent  layer  m  with 
thickness  ALm.  The  wind  has  the  effect  of  shifting  the  turbulent  layer;  therefore,  turbulent 
layer  m  is  moved  by  wind  vector  Vw .  There  are  two  reference  frames  included  in  this 
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figure,  the  first  being  (x,z),  which  is  coincident  with  the  aperture  at  time  ti  and  has  the 

Z-axis  normal  to  the  aperture.  Likewise,  frame  (u,w)  is  coincident  with  the  aperture  at 

time  t2  and  the  W-axis  is  normal  to  the  aperture. 

Figure  7  was  shown  here  simply  to  display  an  uncluttered  drawing.  Figure  8 
shows  the  full  ABL  engagement  geometry  first  introduced  in  [19]. 


Figure  8.  ABL  engagement  geometry. 


This  figure  is  the  same  as  Figure  7  with  all  of  the  necessary  geometrical  constructs 
included.  In  particular,  there  are  three  reference  frames  of  interest: 

1.  Frame  1  =  (ix  Az)  with  origin  Oi  coincident  with  the  aperture  and  the  ix  axis 
aligned  with  the  aperture  at  time  tj. 
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2.  Frame  2  =  (iu ,  iw  j  with  origin  CF  coincident  with  the  aperture  and  the  iu  axis 
aligned  with  the  aperture  at  time  t2. 

3.  Frame  3  =  (iY  ,  L  )  with  origin  Om  located  on  the  mth  turbulent  layer  and  the  ix 
axis  aligned  with  the  turbulent  layer. 

The  notation  i*  is  used  to  denote  a  unit  vector  and  the  subscripts  ‘FI’,  ‘F2\  and  ‘Fm’ 
denote  vectors  written  in  Frame  1 ,  Frame  2,  and  Frame  m,  respectively.  Without  loss  of 
generality,  the  vectors  ix  and  iu  lie  in  the  plane  of  the  aperture  at  times  tj  and  t2  and  VA 

and  VT  are  in  the  (ix,iz)  plane.  In  other  words,  this  is  a  planar  engagement.  The  planar 

engagement  has  been  selected  because  the  ABL  operates  on  missiles  in  boost  phase. 

Hence,  the  aperture  can  be  aligned  so  that  the  majority  of  motion  occurs  in  the  x- 

direction.  Although  this  is  somewhat  restrictive,  it  is  more  than  adequate  for  the 

development  of  the  framework  for  the  application  of  AO  to  the  ABL.  Moreover,  the 

following  notation  is  used  to  denote  the  kinematic  variables  of  the  engagement  geometry: 

y(ti),  y(t2)  =  angles  from  the  ix ,  iu  axes  to  the  turbulent  layer  at  times  ti,  t2- 

d(tj)  =  distance  from  the  aperture  to  the  turbulent  layer  at  time  tj,  measured 
normal  to  the  aperture,  i  =  1,2. 

Lx(tj)  =  distance  from  the  aperture  normal  to  the  tip  of  the  missile,  i  =  1,2. 

Lj(tj)  =  distance  from  the  tip  of  the  missile  to  the  aperture  measured  along  the 
aperture  normal,  i  =  1,2. 

L2(tj)  =  distance  from  the  intersection  of  the  aperture  normal  and  the  missile  to  the 
aperture  measured  along  the  aperture  normal,  i  =  1,2. 

0(tj)  =  angle  measured  from  the  aperture  velocity  vector  to  the  aperture  normal 
vectors,  Z  and  W,  respectively,  i  =  1,2. 
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ALm  =  thickness  of  the  m,h  turbulent  layer. 

Vw  =  turbulent  layer  wind  vector,  yw  =  angle  that  the  turbulent  layer  wind  vector 
makes  with  iv  . 

0(,)(Rx,,t,),  <f>(l)(Ru,  ,t2)  =  inbound,  or  beacon,  wavefront  phase. 
O(0)(Rx2,t,),  <b(0)(Ru2,t2)  =  outbound,  or  high-energy  laser,  wavefront  phase. 
For  the  upcoming  development,  it  will  be  necessary  to  transform  vectors  written  in 
reference  frames  1  and  2  to  vectors  written  in  frame  m.  To  accomplish  this,  two 
transformations  are  necessary.  The  first  transformation  embodies  the  novel  modified 
frozen  flow  hypothesis  [20]  which  takes  into  account  translation  as  well  as  rotation  of  the 
system  as  time  progresses  from  ti  to  t2.  In  mathematical  terms,  given  a  vector  r  ,  written 
in  frame  2,  the  required  transformation  to  express  this  vector  in  frame  1  is 
cosAG  0  sinAG  sin  G(t , ) 

0  1  0  rF2+  0  |vA|At  =  ROT(AG)fF2  +  TRANS  1  (22) 

-sinAG  0  cosAG  cosG(t,) 

where  AG  =  0(ti)  -  Gfe),  At  =  t2  -  tj ,  ROT(AG)  denotes  a  rotation  matrix,  TRANS  1  is  a 

translation  vector  to  frame  1,  and  |«|  denotes  the  magnitude  of  the  corresponding  vector. 

The  reason  for  the  opposite  conventions  in  the  definitions  of  At  and  A0  is  that  a  positive 

At  will  yield  a  positive  A0. 

The  second  transformation  is  as  follows:  given  a  vector  ?  ,  written  in  frame  1,  the 

required  transformation  to  express  this  vector  in  frame  m  is 

cosy(t,)  0  siny(t,)  -d(t,)siny(t|) 

0  1  0  rF1  +  0  =ROT{y(t,)}rF1  +  TRANS  m. (23) 

—  sin y (t j )  0  cosy(tj)  -d(t!)cosy(t,)_ 
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3.3  Aperture  Vector  Projections 

In  this  section,  the  projections  of  the  aperture  vectors,  RXj  and  RGj,  i  =  1,2,  into 
the  turbulent  layer  will  be  determined.  Recall  that  the  only  contributions  to  the 
correlation  between  inbound  and  outbound  wavefronts  are  from  the  layers.  Hence,  the 
spatial  separation  between  aperture  vectors  projected  into  the  layers  must  be  evaluated. 
Figure  9  is  used  to  perform  these  projections.  Figure  9  shows  the  aperture,  target,  and 
turbulent  layer  at  one  time  instant  along  with  the  turbulent  layer  normal  vector,  in .  It  is 
desired  to  determine  an  expression  for  the  point  (X,Y,Z),  in  terms  of  the  kinematic 
variables,  which  requires  evaluation  of  the  line  that  connects  the  points  (x,y,0)  and 
(LT,0,L(t)). 

To  begin,  the  equation  of  the  plane  in  which  the  turbulent  layer  resides  must  be 
determined.  The  normal  vector  to  the  turbulent  layer,  written  in  frame  1  coordinates,  can 
be  expressed  as 

in  =  [-siny  0  cosy]T.  (24) 

Then,  the  equation  of  the  plane  in  which  the  turbulent  layer  lies  is 

-  sin  y(X  -  0)  +  0(  Y  -  0)  +  cos y(Z  -  d)  =  0  =>  -  tany  X  + (Z- d)  =  0 .  (25) 
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Figure  9.  Aperture,  target,  and  turbulent  layer  for  determination  of  projected  vectors. 


Now,  letting  ‘sf’  be  a  scale  factor,  differences  in  the  three  components  for  the  line 
between  (x,y,0)  and  (LT,0,L(t))  are 


X  -  x  =  (sf)(LT  -  x) 

Y  -  y  =  -(sf  )y  (26) 

Z  =  (sf)L(t)  . 


Using  these  expressions  in  the  equation  of  the  plane,  eq.  (25),  yields 


sf  = 


x  tan  y  +  d 

L(t)  -  LT  tan  7  +  x  tan  y 


(27) 


Substituting  eq.  (27)  into  eq.  (26)  gives 


[X  Y  Z]T  = 


Lxd  +  x{L(t)-d} 
L(t)-  tany(LT  -  x) 


L(t)-LT  tany-d 
L(t)-tany(LT  -x) 


L(t)d  +  L(t)x  tan  7 
L(t)-  tany(LT  -  x) 


•  (28) 
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Therefore,  the  projection  of  aperture  point  (x,y,0)  onto  the  mth  layer  is  given  by  (X,Y,Z) 
in  eq.  (28).  In  a  similar  manner,  projections  for  the  vectors  Rx,  and  Rup  i  =  1,  2  can  be 
computed.  Utilizing  the  kinematic  variables  defined  in  Figure  8  yields  the  following 
expressions  for  the  projections  of  the  two  aperture  vectors  at  each  time  instant: 


X  kpFl  — 


npF2 


Xxp 


"kp 


u, 

Vn, 

w„, 


Lk(t,)  — tanvltjtl-TltjSfc,  -Rxk} 


L n (t 2 )  tan y(t2 ){lt (t2 )8 n ,  Ru „ } 


LT(tiWt,)8kl  +{Lk(t,)-d(t,)}Rxk 
{Lk  (t  i )  ~ tan  y(1  i  )LT(t  i  )8k|  -d(t,)}Ryk 
MtiMO  +  Lk(ti)tan'Kti)Rx* 

'  LT(t2)d(t2)8„,+{L„(t2)-d(t2)}Run 
{l„ (t 2 ) — tan y(t2)LT(t2)5nl  -d(t2)}Rvn 
(t2  )d(t2 )  +  Ln  (t2 )  tan  y(t  2  )Ru , 


(29) 


(30) 


where 

~Rxk" 

=  Rxk , 

1 

& 

c 

3 

_ 1 

RYk. 

k  7 

LRVnj 

=  Run ,  the  subscript  ‘p’  denotes  a  projected  vector,  the 

subscripts  ‘FI’  and  ‘F2’  denote  vectors  written  in  frames  1  and  2,  Ski,  8ni  are  Kronecker 

fl  if  k  =  1  fl  if  n  =  l 

deltas  such  that  5kl  =  -J  o  ,c  t  t ,  5nl  =  (  „  ,  and  k  =  1,  2,  n  =  1,  2.  It  can  be 


0  if  k*l 


0  if  n  *  1 


seen  that  eq.  (30)  contains  kinematic  variables  defined  at  time  t2-  Hence,  relationships 
between  the  time  t2  and  time  t|  variables  must  be  determined  so  that  expressions  in 
subsequent  chapters  only  contain  time  ti  quantities. 


3.4  Calculation  of  d(t2),  ft2),  and  the  Apparent  Thickness  of  the  Turbulent  Layer 

In  this  section,  relationships  between  time  t2  kinematic  variables  and  their  time  t| 
counterparts  will  be  determined.  It  should  be  pointed  out  that  d(t2)  and  yfe)  are 
determined  by  the  engagement  geometry  and  the  time  interval  t2  -  tj.  Letting  co  be  the 
slew  rate  of  the  aperture,  it  is  clear  that 
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e(t2)  =  0(t1)  +  ®(t2-tl)- 


(31) 


A  constant  slew  rate  has  been  assumed  since  the  sample  period  is  short. 

To  aid  in  the  calculation  of  y(ta),  Figure  10  will  be  used.  Figure  10  shows  the 
turbulent  layer,  the  aperture  at  time  t2,  the  time  ti  aperture  normal  vector,  Z,  lines  which 
are  normal  to  the  Z-  and  W-axes,  and  some  other  geometrical  constructs. 


TIME  t2 

Figure  10.  Time  t2  geometry  for  calculation  of  YO2). 

As  can  be  seen  in  Figure  10, 
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v(t2)  =  l(t,)  +  p.  (32) 

Utilizing  triangle  QRS  and  the  fact  that  the  sum  of  the  included  angles  is  180  produces 

b  +  90o+e(t,)-e(t2)=180”  =>  b  =  90°  -0(tj)  +  0(t2).  (33) 

Then,  looking  at  the  Z  axis  at  point  R,  the  following  relationship  becomes  apparent: 

b  +  p  +  90°  =  180°  =>  p  =  0(t,)-0(t2).  (34) 

Using  eq.  (34)  in  eq.  (32)  yields  the  desired  result: 

y(t2)  =  Y(t,)+e(t1)-0(t2).  (35) 

In  order  to  determine  dfo),  Figure  1 1  will  be  used.  This  figure  displays  the 
aperture  and  target  at  both  time  instants  along  with  other  necessary  quantities. 
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The  first  step  in  this  process  is  to  determine  the  length  A.  This  can  be  accomplished 
using  the  law  of  sines  on  triangle  BCD  which  produces 


d(ti)  _  A 

sinf^tj  +  ^t,) -90°  sin[90°  - Y(t,) 


d(tj )  cosy(t, ) 

A  —  - ( 

cos[y(  t , )  -f-  e(  t , )] 

Using  the  law  of  sines  once  again,  this  time  on  triangle  BEF,  gives 


d(t2) 

A  + 

vA 

(t2  ll) 

sin 

v(t,)+e(t,)-90-; 

sin[ 

90° 

-yfe) 

Substituting  eq.  (37)  into  eq.  (38)  and  rearranging  yields  the  desired  result: 

_  d(t , )cosy(t1 )  —  |vA |(t2  -tJcos^t^  +  eQ,)] 
cos[e(t2)-9(t1)-Y(t,)] 


(36) 


(37) 


(38) 


(39) 


The  last  quantity  of  interest  in  this  section  is  the  apparent  thickness  of  the 
turbulent  layer.  The  thickness  of  the  layer  is  denoted  by  ALm.  However,  since  the  rays 
from  the  beacon  and  HEL  do  not  pass  perpendicular  to  the  turbulent  layer,  the  apparent 
thickness  of  the  turbulent  layer,  that  is,  that  amount  of  each  slab  which  is  traversed  by 
each  ray,  must  be  determined.  Figure  12  will  be  used  to  evaluate  this  quantity.  This 
figure  shows  the  aperture,  target,  and  turbulent  layer  at  one  time  instant.  Other  quantities 
displayed  in  this  figure  are  the  angles  h  and  f  and  the  apparent  thickness,  denoted  by 
AT  Using  this  figure,  angle  h  can  be  obtained  as 
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(40) 


h  =  tan 


1,(0 

LT(t,)-Rx, 


The  expression  for  h  can  be  simplified  by  considering  the  physical  dimensions  of  the 
problem.  In  particular,  the  ABL  operates  with  Li(t j)  on  the  order  of  100  kilometers, 

Lj(ti)  ~  2  meters,  and  R|x,|  <  0.875  meters.  Therefore,  Lj(ti)  »  Lr(ti)  -  R|x,| .  Hence,  h 


~  n  /  2.  Now,  using  Figure  12,  it  can  be  seen  that  h  =  f  +  y(ti)  =>  f  ~  n  /  2  -  y(tj).  Then, 
by  simple  geometry, 


sinf 


(41) 


Z 


Figure  12.  Engagement  at  one  time  instant  for  apparent  layer  thickness. 
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In  a  similar  manner,  the  apparent  thicknesses  for  the  other  rays  can  be  determined.  The 
expressions  are  the  same  as  that  in  eq.  (41)  with  y(ti)  replaced  by  yfe)  for  time  t2 
thicknesses. 


3.5  Wind  Effect 

To  incorporate  the  effect  of  wind  on  the  turbulent  layer,  the  wind  vector  is 
resolved  into  components  written  in  reference  frame  2.  Thus, 

-  n  -  i  i  — ►  i  l  T 


Vw  = 


v, 


w 


COST)  0 


V, 


w 


sinri 


(42) 


where  T|  =  y(to)  +  yw  and 


V. 


w 


COST) 


sin  rj  are  tne  components 


axis  directions,  respectively.  The  wind  vector  is  written  in  terms  of  frame  2  coordinates 
because  the  effect  of  wind  between  times  t|  and  ti  is  of  interest.  Obviously,  for  t  >  t2, 
wind  effects  also  occur,  but  only  two  time  instants  need  to  be  considered  in  the 
development.  The  results  can  then  be  propagated  to  later  time  instants.  In  the  upcoming 
development,  the  subscript  ‘Vw’  on  a  variable  is  used  to  indicate  that  the  wind  effects 
have  been  incorporated  into  that  quantity. 

First,  consider  the  wind  in  the  W-axis  direction.  Figure  13  shows  the  geometry 
for  this  case.  In  this  figure,  the  aperture  and  target  (source)  are  shown  along  with  the 
turbulent  layer  at  times  ti  and  b-  The  term  source  is  used  since  the  tip  of  the  target,  i.e., 
the  point  of  reflection  of  the  beacon  laser,  can  be  considered  a  point  source.  Notice  that 
the  W-axis  component  of  the  wind  moves  the  turbulent  layer  only  along  the  W-axis.  The 
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two  quantities  that  need  to  be  evaluated  are  the  new  dfo)  distance,  denoted  d(t2  )v  ,  and 


the  new  projected  u  vector,  denoted  uipVwF2 . 


Figure  13.  Wind  effect  in  W-axis  direction. 


It  can  be  seen,  from  Figure  13,  that 


uiPvwF2  =uipF2+Au 


(43) 


where 


Au  = 


0  0 


V, 


w 


sinij 


At 


(44) 


and 


d(t2) 


v  w 


=  d(t2)+  Vw  sinri  At . 


(45) 
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Now,  the  U-axis  contribution  of  the  wind  can  be  examined.  Figure  14  shows  the 
geometry  for  this  case.  In  this  figure,  the  aperture  and  source  are  once  again  shown,  but 
now  the  turbulent  layer  has  moved  from  point  B,  at  time  tj,  to  point  C  at  time  t2.  As  in 
the  case  of  the  W-axis  wind,  the  two  quantities  of  interest  are  d(t2)y  and  uipVwF2 . 


Clearly,  eq.  (43)  holds,  but  eq.  (44)  becomes 


A  u  = 


V, 


w 


COST)  0  0 


rAt 


(46) 


and  the  new  d(t2)  distance  is 


4'A.  =<J(t2)-Ad(t2) 


(47) 


where  Ad(t2)  is  determined  by  simple  geometry  using  triangle  ABC  in  Figure  14: 


Ad(t2)  =  tany(t2)  Vw  cosr|At . 


(48) 


Figure  14.  Wind  effect  in  U-axis  direction. 
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Putting  these  results  together,  that  is,  using  eqs.  (43-48)  leads  to  the  desired  results: 


U,pvw  F2  —  UipF2 


Vw  I  COST)  At 
0 

Vw|sinr|At 


,i=l,2 


(49) 


and 


d(t2)vw  =d(t2)-tany(t2)Vw  cosr|At-i-  Vw  sinr|At 


(50) 


3.6  Projected  Aperture  Vector  Separation  and  Relationships  Between  Aperture 

and  Projected  Aperture  Vectors 

To  begin,  first  consider  the  projected  aperture  vector  separation.  As  previously 
stated,  in  using  a  layered  atmospheric  model,  the  correlation  between  inbound  and 
outbound  wavefronts  is  due  only  to  the  turbulent  layers.  Hence,  the  separation  between 
aperture  vectors,  projected  into  the  slabs,  must  be  determined.  The  notation  used  is  as 
follows:  xkpFm  ,unpVwFm  are  the  projections  of  aperture  vectors  Rxk  ,Run  written  in 

frame  m.  Figure  15  will  be  used  in  this  case.  Here,  the  projected  vectors  are  shown, 
written  in  frame  m  and  frame  2.  From  this  figure,  it  is  clear  that  the  difference  in 
projected  vectors  is  simply 


^npvwFm  ^kpFin  * 

Notice  that  in  order  to  obtain  unpVwF  ,  the  transformations  described  in 
must  be  applied  to  unpVwF2 . 


(51) 

eqs.  (22)  and  (23) 
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Figure  15.  Projected  vector  separation. 


In  the  upcoming  development,  an  expression  for  the  correlation  functions  will 
contain  projected  and  non-projected  aperture  vectors.  In  order  to  determine  the 
correlation  contribution  from  each  layer,  it  is  pertinent  to  have  only  projected  vectors  in 
this  expression.  Hence,  a  relationship  between  projected  and  non-projected  aperture 
vectors  must  be  determined.  First,  consider  the  relationship  between  RxkFI  and  xkpFm. 

Applying  the  transformation  in  eq.  (23)  to  eq.  (29)  gives 


vkpFm 


Mb) 


Lk  (t , )  —  d(t, ) 


cosy(tj) 

0 


Lk(t.)-d(t.) 


RxkH  + 


MhWi) 

cosy(t, ) 

0 


«k. 


=  77~{AxkR*kFl  +Bxk}  >  k  =  1>2 
Lk(tl ) 


(52) 
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Then,  it  is  clear  that 


R*kFl  “^xkl^k^l^kpFm 


(53) 


The  expressions  in  eq.  (52)  were  obtained  using  an  approximation  which  produced  a 
linear  relationship  between  RxkFI  and  xkpFm.  In  the  full  expression,  the  relationship 

between  RxkF1  and  xkpFm  is  nonlinear;  however,  letting 

Lk(t]) »  tany(tl){LT(t,)8k,  —  R|xkFI|| ,  k  =  1,2,  (54) 

the  linear  result  in  eq.  (52)  is  produced.  For  the  case  of  RunF2  and  unpVwFm ,  the  same 
reasoning  has  been  employed  to  produce 
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U 


npvwFm 


Ln(t2) 
cos  y(t2) 

o 


+  a 


Ln(t2)  d(t2)\ 


nF2  + 

_  0  J 

_ 

1 


(f2  ) 


{AunRCnF2+Bun}  ’  n=1’2 


(55) 


where 


a  =  -cosy(t2)d(t2)Vw  -tany(t2){x}. 


(56) 


6un  _  d(t2)vw{COS"f(t2)CT(t2)^nl  +  Sm  y(t2  )Ln  (t,  )}  C„(t2){%},  (57) 


and 


%  =  cosy(t2)Vw  cos ri A t  +  sin y( 1 2 ) Vw  simyA t  —  sin[y(t, )  +  9(t , )]| VA  At  +  d(t , ) sin y(t , ) ,(58) 


n  =  1,2.  Therefore,  it  can  be  seen  that 


R^nF2  ~  ^unj^n  (*2  ^ 


c  -B 

npv  w  Fm  un 


}  ,  n  =  l,2  . 


(59) 


It  should  be  noted  that  Axk,  k  =  1,2,  and  Aun,  n  =  1,  2  are  diagonal  matrices.  Hence, 
multiplication  of  these  matrices  with  each  other  will  be  commutative. 


46 


3.7  Summary 

In  this  chapter,  the  aperture  vectors  have  been  projected  into  the  turbulent  layers, 
the  desired  vector  separation  has  been  determined,  and  the  kinematic  variable 
relationships  have  been  evaluated.  At  this  point,  all  of  the  required  geometrical  analysis 
for  this  problem  has  been  performed.  It  is  now  necessary  to  evaluate  the  correlation  of 
the  Zernike  polynomial  expansion  coefficients  using  the  analysis  in  this  chapter.  As  can 
be  seen  in  Figure  8,  the  inbound  (beacon)  and  outbound  (high-energy  laser)  wavefronts 
propagate  through  different  regions  of  the  atmosphere  at  different  times.  In  other  words, 
the  sensed  wavefront  perturbations  and  the  high-energy  laser  wavefront  are  associated 
with  different  observation  directions.  Therefore,  the  phase  perturbations  associated  with 
the  two  wavefronts  are  different.  This  performance  degradation  is  referred  to  as  an 
anisoplanatic  effect.  To  account  for  anisoplanatic  effects,  the  correlation  between 
wavefronts  arising  from  the  two  directions  must  be  evaluated.  In  Chapter  4,  an 
expression  for  this  correlation  is  derived. 
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IV.  Atmospheric  Modeling 


4.1  Introduction 

The  correlation  between  inbound  and  outbound  wavefronts  is  determined  in  this 
chapter.  Since  measurements  in  the  ABL  scenario  are  of  the  inbound  wavefront’s  phase, 
the  correlation  between  the  inbound  and  outbound  wavefronts  must  be  evaluated  in  order 
to  determine  the  effects  of  turbulence  on  the  outbound  wavefront’s  propagation  path.  The 
outbound  propagation  path  is  important  since  the  HEL  is  fired  through  this  section  of 
atmosphere.  Once  the  correlation  functions  are  determined,  this  data  will  be  used  to 
synthesize  a  stochastic,  linear,  dynamical  system  which  will  serve  as  the  model  for  the 
atmosphere. 

4.2  Correlation  of  Inbound  and  Outbound  Wavefronts 

With  the  analysis  from  Chapter  3  in  place,  the  correlation  between  Zernike 
polynomial  phase  expansion  coefficients  can  now  be  determined.  Much  of  the 
development  of  the  correlation  functions  follows  the  work  of  Takato  and  Yamaguchi 
[26],  Valley  and  Wandzura  [30],  and  Whiteley  [35,  36].  Modifications  have  been  made 
including  the  modified  frozen  flow  hypothesis,  appropriate  layer  thicknesses,  and  the 
ABL-specific  geometry. 

Utilizing  the  Zernike  polynomial  expansion  in  eq.  (17),  the  inbound  and  outbound 
wavefront  phases  can  be  written  as 

*‘”(R*in.t,)  =  Xa<ril(t,)zf(x,F|),  0W(Rx2F,.I,)  =  5>“(t,)Zr(x2FI)  (“t 

f  f 
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and 


0<‘> (RO 1F2 , 1 2 )  =  X af (t 2 )z j (u , F2 ) .  4>w(Ru2F2,t2)  =  XaS°>(t!)Zj(uIP2),  (61) 
i  i 

where  af(l)(t|),  af(0)(ti),  aj<l)(t2),  and  aj(0)(t2)  are  the  Zernike  polynomial  phase  expansion 
coefficients.  The  superscripts  (i)  and  (o)  refer  to  the  inbound  and  outbound  wavefronts. 
The  expansions  used  in  this  work  will  include  two  modes,  that  is,  f  =  2,  3  and  j  =  2,  3. 
Zernike  mode  number  one  is  termed  piston  (aperture  averaged  phase)  and  is  not  a 
distortive  contributor.  Hence,  it  is  ignored  in  the  expansions.  The  general  forms  in  eqs. 
(60)  and  (61)  will  be  used  to  allow  the  framework  to  be  built  using  as  many  modes  as  one 
desires.  The  decision  to  use  two  modes  will  be  discussed  later  in  this  work. 

These  equations  can  be  written  in  more  compact  forms  as 

4>(l"(RxH,„t,)  =  £a?*>(t,)Zr(iim),k  =  l,2  (62) 

f 

and 

O(*'0)  (Ru nF2  5 1 2 )  “  a|i,0) (t |  )Zj (u np2 ) ,  n  =  1,2  (63) 

j 

where  the  notation  0(1,0)(*,«)  and  a(l,0>(t)  designate  the  inbound  and  outbound  phase  and 
Zernike  coefficients  depending  on  which  aperture  vectors  are  used.  In  other  words,  if 
using  aperture  location  x1F1  in  eq.  (62),  i.e.,  k  =  1,  then  the  inbound  wavefront, 
superscript  (i),  is  under  consideration.  If  aperture  location  x2FI  in  eq.  (62)  is  used,  i.e.,  k 
=  2,  then  the  outbound  wavefront,  superscript  (o),  is  under  consideration.  The  same  holds 
for  the  expressions  in  eq.  (63).  This  notation  allows  derivation  of  the  correlation  between 
inbound  and  outbound  wavefronts  to  be  evaluated  in  the  general  case,  with  the  desired 
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combination  of  wavefronts  substituted  later.  The  Zernike  coefficients  can  now  be 
determined  as  specified  in  eq.  (18): 

alri'o,(t1)  =  jw(xkF1)zf(xkF,)4>0»,(RxkFI,t,)dxlFI  (64) 

and 

af°’(t,)  =  |w(uJF2)zj(u-R)<t.«'">(RunF;,t2)dii„R  (65) 

where  the  integration  is  taken  over  the  entire  two-dimensional  plane  and  W(»)  was 
defined  in  eq.  (19).  Since  it  is  desired  to  determine  the  correlation  between  the 
wavefronts,  and  the  wavefronts  have  been  expanded  using  eqs.  (60)  and  (61),  it  follows 
that  the  desired  expression  can  be  evaluated  by  computing  the  correlation  between  the 
Zernike  expansion  coefficients.  Hence,  the  correlation  becomes 

r.rc., =E{a“'0l<t,)afo>(12)} 

=  E{JJw(xlF,)zf(xkFI)w(u,F2)zr(s,F2)o<l">(Rxkp,,tl)®<l’")(Ru„P2,t2)dxkF,du„F2} 

=  jjw(xkFI)z,(xkF,)w(unF2)Zf(u,F2)E{<I)<l'",(RxlF,,t,)<l>,i”>(RunF2,t2)}dxkFldu„F2 


(66) 

where  the  last  expression  was  obtained  using  the  linearity  property  of  the  expectation 

operator  and  T  (i0),  ,  (l0),  ,  is  the  correlation  between  Zernike  coefficients. 

r  ap •  >  ( t , ) a *  '  ( 1 2 ) 


In  order  to  evaluate  the  integral  in  eq.  (66),  it  is  necessary  to  determine 

E{on"l(RxkFI,tl)d>"«(RiinR,t,)}.  (67) 


Some  of  the  properties  described  in  Chapter  2  are  now  necessary.  In  particular,  using  a 
layered  atmospheric  model  and  assuming  that  all  layers  are  of  equal  thickness  and  the 
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index  of  refraction  fluctuations  are  uncorrelated  layer  to  layer  (by  suitably  placing  the 
layers  as  discussed  in  Chapter  2),  a  wavefront’s  turbulence-induced  phase  can  be 
expressed  as 

0(x)  =  KSALm>nm(^t)  (68) 

m 

where  k  =  2k  /  X  is  the  spatial  wavenumber,  X  is  the  wavelength,  m  is  a  layer  subscript, 

AT  ^  is  the  apparent  thickness  of  the  mlh  layer  as  seen  by  the  kth  wavefront,  and  nm(x,t) 

is  the  mlh  layer’s  randomly  fluctuating  index  of  refraction  term.  Using  eq.  (68)  in  eq.  (67) 
produces 

E|<J)(i'0)  (Rx  k  F1  )o(i,0)  (Ru  n  F2  )|  =  e|k  2  ^  AL(k)  AL(^  n  m  (x  k  pF1 ,  t ,  )n  m  (u  n  pF2 , 1 2  )|  (69) 

where  ALm(k),  ALm(n)  are  the  apparent  layer  thicknesses  for  the  kth,  nth  wavefronts  and 
x k pFj  ,unpF2  are  the  projected  aperture  vectors  written  in  frame  1  and  frame  2.  Frame  2  is 

used  for  unp  since  the  modified  frozen  flow  hypothesis  (FFH)  has  not  been  included.  In 
fact,  applying  the  modified  FFH,  as  described  in  eq.  (22),  yields 

nm(^npF2 ’U)  =  nm  (QnpFl ’^1  )  •  (^0) 

In  other  words,  the  modified  FFH  takes  into  account  target  and  source  motion  and  thus 
relates  frame  2,  time  t2  variables  to  frame  1,  time  ti  variables. 

Now,  transforming  the  vectors  in  eq.  (69)  to  frame  m  and  incorporating  the  effects 
of  wind  produces 
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(71) 


E|o(i,0)(RxkF1)<I)<i’0)(RunF2)|  =  E|k2^  AL(1^)AL(1^)nm(xkpFm,tl)nm(unpV(tFm,t1)| 

^[K2AL2mE{nm(xkpFm,t,)nm(unpVwFm,t,)} 

_ __m _ _ 

cosy(t,)cosy(t2) 

where  the  cosine  terms  arise  because  the  expression  for  the  apparent  layer  thickness  has 
been  included.  It  can  be  seen  that  the  term  in  brackets  is  the  spatial  autocorrelation  of  the 


phase  in  the  mth  layer.  Thus 


E{<I>'‘°’(RxliF1)o‘i'“'(RunF,)}  = 


2 (*  k  pFm  ’  ^  n  pvwFm  ’  *  1  ) 
m _ 

cosy(t>)cosy(t2) 


(72) 


where  t , )  denotes  the  spatial  autocorrelation.  Under  the  assumptions  of 

homogeneous  and  isotropic  statistics,  as  addressed  in  Chapter  2,  ro  ( t , )  will  be  a 
function  only  of  the  vector  separation  between  its  arguments.  Therefore,  using  eq.  (51), 

(*kpFm  >^npvwFm  »  C  )  —  ’  C  )  —  npv„,Fm  *  kpFm  4  ]  )  •  (23) 

Substituting  eq.  (73)  into  eq.  (72),  using  this  result  in  eq.  (66),  and  interchanging  the 
order  of  integration  and  summation  produces 

Vo,*!' . „„=E{a“«,)a“(t2)} 


-I 


cosy(t,)cosy(t2) 


(^npvwFm  *k pFm  )w(xk FI  )^f  (* k  FI  )^(u  „  F2  )Zf  (u nF2  )dXkF1du nF2  • 


(74) 
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Due  to  the  fact  that  eq.  (74)  contains  projected  and  non-projected  aperture  vectors,  the 
transformations  given  in  eqs.  (53)  and  (59)  must  be  used  to  obtain  only  projected  vectors. 
Applying  these  transformations  yields 


. (l,  )a'J' "’02)  X 


1 


cy(t,)cy(t2) 


Jjr0m(5npv..Fm-ApFJW 


kpFm 


■M 


*w 


npvwFm 
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A"1  r 

ML(t2)u 
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npvwFm 


Bd 


A~'L(t, ) 
R 


ldXkPFm 


Au'L(t2) 


R 


du 


npvwFm 


(75) 


where  Ax'  ,L(t,),Bx  are  Axl'  ,Ll(t,),Bxl  if  k  =  1  or  Ax2  ,L2(t,),Bx2  if  k  =  2,  and 
likewise  for  the  u  components. 

There  are  well-known  forms  for  the  Power  Spectral  Density  (PSD)  corresponding 
to  the  spatial  autocorrelation  of  the  phase,  namely,  the  Kolmogorov  and  von  Karman 
spectra.  Since  the  autocorrelation  of  the  phase  and  its  PSD  form  a  Fourier  transform  pair, 
transform  techniques  can  be  applied  to  eq.  (75)  to  obtain  a  PSD  term.  This  analysis 
follows  the  development  in  the  work  of  Takato  and  Yamaguchi  [26]  and  simplifies  the 
integrations  described  in  eq.  (75).  The  following  transform  property,  among  others,  will 
be  useful,  assuming  that  s(r)  is  Fourier  transformable:  FT{*}  denotes  the  Fourier 
transform  of  { • } . 

•  Scaling  and  shifting  property  [7]: 


FT{s(a,x  +  b,y  +  c,  ,a2x  +  b2y  +  c2)j 


-|2ic(xu<!),+y0m2 


a7  -b,  a, 

— l($t, — -co.  +  — CD 

D  2  D  D 


(76) 
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,  _  ,  ,  b,c, -b9c,  ,  a,c, -a,c7 

where  D  =  a,b2  ~a2b,  ,x0  =  —  —  - — ,  and  y0  = 


D  D 

First,  consider  the  integral  over  u  Fm  in  eq.  (75).  Letting 


A..  = 


AuU  0 
0  A 


u2,2 


.»„  = 


B 


ul.l 

0 


,  u 


npvwFm 


UnpvwFm(^) 
^  npvwFm 


,  k  = 


(77) 


the  terms  inside  the  weighting  function  and  Zernike  polynomial  can  be  expressed  as 

A..1  „  1  r  f  k>(f  2  )Unpv„,Fm  (^)  Bu]  , 


R 


|L(t2)unpv  Fm  Buj  — 


^ °  A  L(t,)n  ,<2) 


-|T 


V  RAuu 


RA 


ul,1  ) 


L2  7wnpvwFmv 
^^u2,2 


(78) 


Using  eqs.  (75)  and  (77),  the  components  ai,  bi,  Ci,  a2,  b2,  c2,  D,  x0,  and  y0  can  be 
identified  and  the  Fourier  transform  of  Zj(®)W(*)  can  be  written  as 


FU  Z: 


Y{Ut2)unpv>Fm-Bu}  W  ^-{L(t,)unpv  Fni  ~B„} 
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L2(t2)  |det(A”') 
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(79) 


where  Qj*(f)  is  the  Fourier  transform  of  Zj  (f)W(f),  ‘  ’  denotes  complex  conjugation,  and 


-j27tic* 


the  dot  (•)  in  e  L<‘2>  denotes  dot  product.  Using  this  result  and  Parseval’s  Theorem 
[8]  in  eq.  (75)  produces 

^a<fi’0)(t1)a(ji-0)(t2)  =  X  cy(t ,  )cy(t  2  )  ^npvwFm  _^kpFm)}W  R  {^(1,  )XkpFm  “BXJ 
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Y{L(t,)xkpFm-Bx} 
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L(t2)  |det(A*‘)| 
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LL(t2)  J 

L  R  J 

ldxkpFmA-'dk. 


(80) 


Using  the  Wiener- Khinchin  theorem  [7]  and  the  shifting  property  of  the  Fourier  transform 
gives 
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FT{r«,.  (u„„.Fn,  -  xlpFm)}  =  W..  (K)e-1“"l-s'"-> 


(81) 


where  W4  (k)  is  the  PSD  of  the  spatial  autocorrelation  of  the  phase  in  the  mth  layer. 
Substituting  eq.  (81)  into  eq.  (80)  and  rearranging  yields 

1  R  1 


cY(t,)cy(t2)  L(t 2 )  |det(A“')| 


jw0m(K)e"j2^L(t2)^Q; 


RAU  _ 

- — K 

L(t2) 


Jej2-Vn,w  ^-{L(t|  )xkpFm  -  Bx}  Zf  ^-{L(t,)xkpFm-Bx} 
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A  L(t , ) 
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|dxkpFm 


A..'dk  -(82) 


Using  the  Fourier  transform  properties  once  again,  the  term  in  parentheses  in  eq.  (82)  can 
be  transformed  to  yield 
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R2 


cy(t,)cy(t2)  L(t,)L(t2) 


|det(Ax)||det(Au)| 


,;jw,n(Ue 


-j27tK« 
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(83) 


where  Ax  and  Au  are  2-by-2  diagonal  matrices.  Thus,  multiplication  of  these  two 
matrices  is  commutative. 

Equation  (83)  provides  the  desired  result,  that  is,  the  correlation  between  Zernike 
polynomial  phase  expansion  coefficients.  However,  double  integrations  are  not  feasible 
in  an  on-line  situation.  Therefore,  it  would  be  desirable  to  obtain  a  simplified  expression 
by  integrating  out  one  component  of  the  frequency  vector.  To  simplify  eq.  (83),  consider 
converting  the  Cartesian  integral  into  an  equivalent  polar  integral.  Letting 

Ax'A“'K  =  k  =>  A^A^'dk  =  die,  (84) 

it  can  be  seen  that  (recalling  the  Ax  and  Au  are  diagonal,  so  that  they  commute) 
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k  =  AxAu*  ,  =  ,  and  ^  ^ 


L(t2 )  L(t2) 


L(t , )  L(t , ) 


(85) 


-*  B  B 

Also,  W*  (ic) ->  W0m(AxAuK)  and  letting  a  =  — -JL-+— yields 


-j27cic»a  _ ^  ^-j2TtAxAuK*d 


Therefore,  eq.  (83)  becomes 


r  =  y _ !. 

a W (*2 )  *  cy( t  tc 


R2 


cy(t,)cy(t2)  L(t,)L(t2) 


|det(Ax)||det(Au)| 
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(86) 


For  notational  simplicity,  let 


ic,  =  ic  =  A  A„k  ,  ic. 


RAA 


— ic  ,  and  ic3  = 


ra:a.,  - 

-K, 


L(t2)  L(t, ) 

where  ic,  =  vector  whose  angle  to  ix  is  0,,  ic2=  vector  whose  angle  to  ix  is  02,  ic3  = 
vector  whose  angle  to  ix  is  03,  and  a  =  vector  whose  angle  to  ix  is  0a.  Thus, 
converting  eq.  (86)  to  an  equivalent  polar  integral  produces 
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Now,  letting 
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it  can  be  seen  that 
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and 
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(91) 


Using  Noll’s  expressions  for  the  Fourier  transform  of  the  Zernike  polynomials  [18],  that 
is, 

Qf(K,0)  =  V^T(-D2V2(l  S""^  n,+l^ — -cos|mf0  +  -^-(l  —  Sm  ,0)[(— 0f  —  l]|  (92) 

where  n,  m  are  the  radial  degree  and  azimuthal  frequency  of  the  Zernike  polynomials, 

J  ,(•)  is  a  Bessel  function  of  the  first  kind  of  order  (nf  +  1),  and  5m  0  is  a  Kronecker 


delta  such  that  8  0 


1  if  mf  =  0 
0  otherwise 


,  in  eq.  (88)  yields 
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Unfortunately,  there  is  not  a  readily  available  closed-form  solution  to  the  integral 
over  0 .  The  closed  form  is  essential  for  the  real-time  applicability  of  this  work.  Without 
performing  the  integration  over  0 ,  numerical  double  integrations  would  be  required  in 
real  time.  Since  the  turbulence-induced  wavefront  deformations  change  relatively  quickly 
and  the  system  must  be  updated  at  speeds  commensurate  with  these  changes,  on-line 
numerical  double  integrations  are  not  feasible.  However,  single  integrations  are 
acceptable  and  hence  it  is  desirable  to  integrate  out  the  0  dependence.  Therefore,  to  aid 
in  the  evaluation  of  the  integral  over  0,  two  additional  assumptions  are  now  introduced. 
The  two  assumptions  are  as  follows: 

1.  y(t,)  =  0°  and  2.  A  t « 1  sec.  (94) 

Recall  that  y(ti)  is  the  angle  from  the  plane  of  the  aperture  to  the  turbulent  layer  at  time  ti. 
The  first  assumption  can  be  justified  by  noticing  that  the  layered  turbulence  model 
requires  y(t|)  ~  0  .  This  would  correspond  to  the  classical  situation  of  an  optical  imaging 
system.  Clearly,  y(ti)  ~  90  is  not  a  practical  physical  situation  since  the  propagation 
paths  would  always  lie  within  a  single  layer.  The  second  assumption  comes  from  the 
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prediction  of  correlation  times  on  the  order  of  tenths  of  milliseconds.  With  these 
assumptions,  it  can  be  seen,  using  eqs.  (31)  and  (35)  that 

0(t,)»0(t2)=»Y(t2)*Y(t,)=>Y(t2)«O-.  (95) 

Also,  using  eqs.  (53)-(59),  (89),  and  (94),  it  is  found  that 

AxU  =Ax2,2  ar,d  Aul ,  —  Au22.  (96) 

Therefore,  after  some  manipulation,  the  integral  over  0  in  eq.  (93)  reduces  to 

INTg  =  —  J  e  j2T'HL0S(e  e“)  cos|(mf +mj)0  +  gj  +  cos|(mf -mjjG  +  hjjdB  (97) 
2  9n 

where 

*  =  {f  {0  -K »!(-')' -■]+('  -s^oXc-D1  - 1]}}  (98) 

and 

h  =  {f  {(1-8™, „)[(-!)'  -  l]-(l  „)[(-l)J  - 1]}}  ■  (99) 

By  a  change  of  variables,  a  trigonometric  identity,  and  the  definition  for  a  Bessel  function 
of  the  first  kind  [9],  denoted  by  J»(*)>  it  is  found  that 

3(m‘+mi)  ,  _  V 

INT§  =  cos|(mf  +  mj)ea  +g}jc(-l)  2  J(%+Bj)(2ick1  |a|) 

+cos{(mf -mj)0a  +h]jt(-l)  2  (2k k,  |d|) .  (100) 

The  details  of  integrating  eq.  (97)  are  shown  in  Appendix  A.  Therefore,  the  expression 
for  the  correlation  between  Zernike  polynomial  phase  expansion  coefficients  becomes 
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3(m'~mi)  , 

+  cos{(mf  -mj)9a  +h|(-l)  2  J[mf_mj|( 


The  final  component  of  eq.  (101)  that  must  be  specified  is  the  power  spectrum.  In  this 
work,  the  von  Karman  spectrum,  eq.  (13),  is  used  and  with  a  weighted  turbulence- 
dependent  layer  model  can  be  expressed  as 

/83 


V3r 


w*  (k,Cm)  =  - 


2(2  7t)3 

Using  eq.  (102)  in  eq.  (101)  produces  the  final  result: 


( 6.883  y/  2  n" 
U-9 \)  1  ' 
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(103) 


4.3  Calculations 

Now  that  the  analysis  has  been  set  into  place,  the  correlation  functions  for  the  first 
few  combinations  of  Zemike  polynomials  will  be  determined.  These  calculations  are 
performed  by  numerically  integrating  eq.  (103)  over  the  spatial  frequency  vector. 
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The  following  scenario,  similar  to  that  used  by  Whiteley,  et.  al.  [35,  36]  is 
considered:  one  turbulent  layer  is  used  and  is  placed  a  distance  of  25,000  meters  from  the 
aperture  at  time  ti.  The  geometry  is  shown  in  Figure  16  and  the  engagement  parameters 
are  as  follows: 

Va=[300  0  o]T ,  VT  =  [o  o  o]T,vw=[o  o  o]T  (104 

L2  (t , )  =  50,000 ,  Lt  (t , )  =  1 , ■ y(t , )  =  0° ,  R  =  0.875 

where  the  velocity  and  distance  units  are  meters/sec  and  meters,  respectively. 


vT=[o,o,o] 


Figure  16.  Engagement  geometry  for  calculations. 


For  the  simulations  and  modeling  to  follow,  two  Zemike  modes  (tilts)  are  used  to 
approximate  the  wavefront  phases.  This  number  may  seem  quite  small,  however,  the 
order  of  the  linear,  stochastic,  dynamical  system  modeling  the  atmosphere  (plant)  that  is 
derived  from  these  correlations  is  two  times  the  number  of  Zernike  modes  used  for  each 
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wavefront  expansion.  Thus,  if  more  Zemike  modes  were  used,  the  order  of  the  system 
would  increase  to  the  point  that  advanced  control  concepts  would  be  difficult  to  apply  in 
an  on-line  application  as  required  in  adaptive  optics.  The  Zemike  modes  used  in  this  case 
are  Z2(x)  and  Z3(x) .  These  modes  represent  wavefront  tilt  at  the  aperture.  Tilt  can  be 
described  as  the  least  squares  sense  best  fit  to  the  aperture  phase.  It  is  an  important 
component  of  wavefront  expansions  since  x-  and  y-tilts  represent  about  87%  of  the  power 
in  the  phase  fluctuations  [23].  The  first  Zernike  mode,  Z,(x) ,  is  termed  piston  and 
represents  the  aperture  averaged  wavefront  phase.  Piston  is  not  a  distortive  contributor 
and  is  not  considered  in  any  calculations. 

The  following  notation,  as  introduced  in  Chapter  3,  has  been  adopted  to  represent 
the  Zernike  polynomial  phase  expansion  coefficients:  an(1),  an(0)  are  Zernike  coefficient  n, 
(n  =  2,  3),  of  the  inbound  (i)  and  outbound  (o)  wavefronts.  Using  two  Zernike  modes, 
there  are  ten  correlation  functions  of  interest.  To  see  this,  consider  the  stacked  vector 
which  contains  the  Zemike  expansion  coefficients: 

a“(t)  =  [a?'(t)  af(t)  «“(t)  a<«(t)]T  (105) 


where  the  superscript  ‘a’  implies  atmosphere.  Integration  of  eq.  (103)  will  yield  the 
correlation  kernel  function  of  the  vector  in  eq.  (105)  which  can  be  expressed  as 


E{a(a)(t)a(a)(t  +  x)T}  = 


"Efafa^}  E{a^°}  E{a<0a!,o)}  E{a^a‘o)}' 

Efa^af}  Efa^a*0 }  Efa^a^}  E{a<°a<0)} 


Efa^an  Ejaj  a2°  )  E{a)0,a^}  E{a^a^} 

Efa^a^}  Efa^a'0* }  E{a<0)a‘0)}  E{a<0)a<o)}_ 


(106) 


where  t  =  At  =  tj+i  -  tj.  This  matrix  is  symmetric,  thus  only  the  diagonal  and  upper  or 
lower  diagonal  terms  need  be  calculated.  Hence,  there  are  ten  correlation  functions  to  be 
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ans 


E{a3i  :  a3i} 


Figure  21.  Eja^a^  J  and  Eja^a^j  versus  At. 


As  can  be  seen  from  these  figures,  four  of  the  correlations  are  at  least  fifteen  orders  of 
magnitude  less  than  the  others.  Therefore,  Eja^a^j,  Eja^a^J,  Eja^a^  j,  and 

E|a20)a(,0>  j  will  be  set  to  zero  in  the  upcoming  work.  Now  that  the  calculations  are  in 

place,  it  is  necessary  to  realize  an  underlying  linear,  stochastic,  dynamical  system  which 
has  correlations  that  are  representative  of  those  in  Figure  17-21.  The  next  section 
addresses  this  task. 

4.4  Stochastic  Modeling 

Since  the  correlation  functions  calculated  in  Section  4.3  are  due  to  the  turbulent 
atmosphere,  the  linear,  stochastic,  dynamical  system  to  be  derived  will  in  turn  be  a  model 
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of  the  atmosphere,  i.e.,  the  plant.  Due  to  the  random  nature  of  the  atmosphere,  this  model 
will  be  stochastic.  Let  the  linear  system  to  be  identified  be  represented  by  [16] 

da(a)(t)  =  Faa(a)(t)dt  +  GadBa(t)  (107) 

where  Ba(t)  is  Brownian  motion  with  statistics  E{Ba(t)}  =  0  and 

E{[Ba(t)-Ba(tO][Ba(t)-Ba(t')]T}  =  jQa(t)dx=JldT  (108) 

t  t 


for  all  t  >  t ,  I  is  the  identity  matrix,  and  Fa,  Ga  have  the  following  forms: 

ffn  f  12  f.3  fi4l  f§i .  §>2! 


f  f  f  f 

l2 1  1 22  l23  l24 

f  f  f  f 

1 31  *32  1 33  1 34 

f  f  f  f 

Ml  l42  1 43  *44 


§11 

§12 

§21 

§22 

§3! 

§32 

§41 

§42 

(109) 


The  dimensionality  of  Fa  was  determined  by  the  number  of  states  in  the  state  vector,  i.e., 
a<a)(t)  =  [a^t)  a3°(t)  a^t)  a30> (t)] .  The  four  element  state  vector  contains  the  x- 


and  y-tilt  Zernike  expansion  coefficients  of  the  inbound  and  outbound  wavefronts.  Tilts 
were  used  since  they  contain  most  of  the  power  in  the  wavefront  aberrations.  Higher 
order  modes  can  be  used,  but  their  inclusion  does  not  significantly  improve  the  developed 
model.  Ga  was  selected  as  a  4-by-2  matrix  through  iterative  testing.  It  was  desired  that 
the  realized  system’s  correlation  kernel  functions  accurately  represent  those  which  were 
calculated.  In  this  case,  a  4-by-2  input  noise  matrix  sufficed.  Notice  that  in  eq.  (108),  the 
diffusion  of  Ba(t)  was  set  to  I.  This  simply  eliminates  one  of  the  design  parameters  and 
allows  system  identification  to  be  performed  with  more  ease.  Also,  it  can  be  seen  that  eq. 
(107)  lacks  a  Bu(t)dt  term.  This  is  due  to  the  obvious  fact  that  the  atmosphere  is 
uncontrollable. 
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Using  eq.  (107),  the  correlation  kernel  matrix  of  the  state  process,  denoted  Pa(t,x), 
satisfies 

E{a(a)  (t)a(a)  (t  +  x)T}  =  Pa  (t,  t  +  x)  (110) 

where 

Pa(t,t  +  x)  =  Pa(t)eF‘\x>0  (111) 

and  Pa(t)  is  the  correlation  matrix.  Taking  the  derivative  of  eq.  (1 1 1)  with  respect  to  x 
and  manipulating  yields 

^4^  =  Pa(t)eFjTFaT  =Pa(U  +  x)FaT.  (112) 

dx 

Equation  (1 12)  is  needed  to  obtain  the  system  matrix.  The  calculations  in  Section  4.3 
provide  the  correlation  kernel  matrix.  The  derivative  of  this  matrix,  with  respect  to  x,  can 
then  be  obtained.  Hence,  the  only  unknown  in  eq.  (1 12)  is  the  system  matrix.  Therefore, 
using  the  correlation  data  in  Section  4.3  and  eq.  (112),  the  system  matrix,  Fa,  can  be 
identified.  In  order  to  obtain  the  input  matrix,  Ga,  consider  the  correlation  of  the  state 
process  which  satisfies  the  following  differential  equation: 

P.(l)  =  F.Pa(t)  +  P.(t)F.T+G.IGl,  (113) 

Letting  t  — >  °°  while  maintaining  a  fixed  separation  between  t  and  t  +  x,  which  is  the 


steady-state  case  and  produces  a  stationary  process,  provides 


where  the  subscript  ‘ss’  implies  steady-state.  Pa  (t)  is  found  from  the  calculated 

correlation  kernel  data  by  letting  x  -  0.  Then,  Ga  can  be  determined  by  taking  the 
Cholesky  decomposition  [16]  of 

G,Gl=-(F,P,Jt)  +  P.Jt)r:).  (115) 

After  performing  the  computations,  the  atmospheric  model  becomes 

-410.35  0  420  0  1  [.4445  0 

0  -437.47  0  464.41  ,  ,  0  .2687 

a(!°(t)dt  + 

-207.88  0  199.75  0  .27  0 

0  -186.53  0  1 65.59 J  0  .2516 

=  Faa(a,(t)  +  GadBa(t).  (116) 

By  computing  the  eigenvalues  of  the  Fa  matrix,  namely,  -29.48,  -70.39,  -181.12,  and 

-201 .48,  it  can  be  seen  that  the  system  is  stable.  Figures  22-24  show  the  nonzero 
correlation  functions  from  the  calculations  in  Section  4.3  with  those  from  the  atmospheric 
model  in  eq.  (116)  superimposed.  These  three  figures  convey  the  fact  that  the  model  does 
provide  a  good  representation  of  the  calculated  correlation  functions.  However,  there  is 
some  modeling  error.  It  would  appear  that  using  more  Zemike  modes  in  the  wavefront 
expansions  would  provide  a  better  representation  of  the  calculated  correlation  functions. 

In  fact,  using  5  Zernike  modes  in  the  wavefront  expansions  was  investigated.  Although 
the  accuracy  did  improve  with  more  modes,  the  modeling  error  was  consistently  present. 
There  is  a  trade-off  between  modeling  accuracy  and  computation  time.  More  modes 
increases  accuracy,  at  the  expense  of  increased  computational  delays.  Since  the  modeling 
error  decreased  only  slightly  with  more  modes,  the  decision  was  made  to  use  tilts  (two 
modes)  in  the  wavefront  expansions. 
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Figure  24.  E{a2o)a20)}  and  Eja^V^j ,  calculated  and  model  versus  At. 


One  final  point  concerning  the  Zernike  phase  expansion  coefficients.  According 
to  Noll  [18],  these  coefficients  are  well  modeled  as  zero-mean,  Gaussian  random 
variables.  They  are  zero-mean  because  the  wavefront  phases,  <1>(1>  (Rx , ,  t , ) , 

<J>(,)(Ru,  ,t2),  <f>(0)(Rx2,t,),  and  <E>(0)(Ru2,t2)  are  the  deviation  of  the  phase  from  the 

ideal,  or  unperturbed  wavefront.  Gaussianness  is  the  result  of  the  wavefront  phases  being 
the  sum  of  a  large  number  of  random  variables.  The  Central  Limit  Theorem  [3]  then 
states  that  the  phases  will  be  Gaussian. 


4.5  Summary 

At  this  point,  the  correlation  between  Zernike  polynomial  expansion  coefficients 
has  been  determined.  Using  these  correlations,  a  state-space  model  of  the  atmosphere 
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was  constructed.  In  order  to  characterize  the  AO  control  system  fully,  a  model  of  the 
mirror  must  be  developed  and  an  output  equation  must  be  specified.  The  next  chapter 
addresses  the  mirror  and  complete  models. 
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V.  Deformable  Mirror  and  Complete  Models 


5.1  Introduction 

In  the  upcoming  work,  a  model  for  the  deformable  mirror  is  derived.  The 
deformable  mirror  (DM)  currently  used  has  a  continuous  mirror  surface  worked  by  49 
actuators.  It  will  be  assumed  that  all  actuators  can  be  driven  by  control  voltages,  i.e., 
there  are  no  slave  actuators.  Actuator  linearity  is  also  assumed  which  implies  that  the 
graph  of  actuator  displacement  versus  applied  control  voltage  is  a  straight  line. 
Additionally,  the  deformable  mirror  is  a  linear  device.  In  this  case,  inter-actuator 
superposition  holds.  For  example,  if  a  voltage  is  applied  to  actuator  1  and  the  mirror’s 
response  is  measured,  then  a  control  voltage  is  applied  to  actuator  2,  the  sum  of  the  two 
responses  would  be  equivalent  to  applying  the  voltages  simultaneously.  After  completing 
the  mirror  modeling,  the  complete  system  will  be  formed  by  appending  the  mirror  and 
atmospheric  models.  Output  equations  for  the  system,  that  is,  measurements  of  the 
inbound  and  reflected  wavefronts,  will  be  covered.  Lastly,  the  continuous-time  system 
will  be  converted  to  an  equivalent  discrete-time  system  for  implementation  in  a  digital 
computer. 

5.2  Steady-State  Mirror  Behavior 

In  this  section,  it  is  desired  to  determine  the  relation  between  applied  control 
voltages  and  the  mirror’s  behavior.  The  term  steady- state  implies  that  a  finite  time  has 
elapsed  since  the  last  occurrence  of  a  control  voltage,  i.e.,  the  mirror  is  not  deforming. 
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The  analysis  used  to  determine  the  mirror  model  and  output  equations  follows  the 
derivation  in  the  work  of  Von  Bokem,  et.  al.  [32]  and  Von  Bokern  [33].  However,  the 
determination  of  some  of  the  required  matrices,  as  detailed  in  Appendix  B,  are  specific  to 
this  work.  This  is  due  to  the  functions  chosen  to  represent  certain  parameters,  as  will  be 
shown. 

To  begin  the  derivation,  consider  the  HEL  beam.  The  HEL  wavefront,  prior  to 
reaching  the  DM,  is  planar  and  in  phase.  For  a  section  of  the  DM  that  has  not  been 
actuated,  the  change  in  Optical  Path  Lengths  (OPL)  is  zero.  However,  there  is  a  decrease 
in  the  OPL  for  a  section  of  the  DM  that  has  moved.  This  decreased  OPL  causes  the 
reflected  light’s  path  length  to  be  decreased  by  two  times  the  change  in  OPL.  Thus,  the 
absolute  phase  of  the  reflected  wavefront,  corresponding  to  a  section  of  the  DM  that  has 
been  actuated,  is  larger  than  for  a  section  of  the  DM  that  has  not  moved.  Since  the 
incident  HEL  wavefront  is  planar,  any  phase  distortion  in  the  reflected  image  will  be 
solely  due  to  the  DM’s  shape. 

To  determine  the  steady-state  behavior  of  the  mirror,  influence  functions  will  be 
used.  Influence  functions  represent  the  effect  of  a  single  actuator  voltage  on  the  mirror 
shape  and  manifest  themselves  in  differences  in  OPL  for  the  HEL.  Utilizing  the  influence 
functions  and  the  fact  that  there  are  in  general  nA  actuators,  the  DM’s  shape  can  be 
represented  by 

0(m)(x)  =  ]>A(x)vak  (117) 

k-1 
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where  x  is  a  spatial  location  on  the  DM,  Ik  (x)  is  the  kth  influence  function,  and  vak  is  the 
control  voltage  applied  to  the  kth  actuator.  Now,  the  reflected  wavefront’s  phase 
distortion,  excluding  piston,  is  given  by 


®'">(5)  =  ^Ik(x)v,k-E  £lk(x)v,k 


k=l 


Lk=l 


(118) 


where  the  expected  value  term,  E[«],  is  piston  and  the  superscript  ‘m’  implies  the  mirror. 
The  reflected  wavefront  can  also  be  written  as  a  linear  combination  of  Zernike 


polynomials,  that  is, 


<b(m)  (x,  t)  =  X  a[m)  (t)Zj  (x) 


(119) 


where  ai<m)(t)  are  the  Zernike  polynomial  phase  expansion  coefficients  associated  with  the 
mirror.  As  shown  in  Chapter  2,  the  mirror  expansion  coefficients  are  given  by 

a[m) (t)  =  JW(x)Zi(x)0<m)(x,t)dx  .  (120) 

Substituting  the  reflected  wavefront’s  phase  distortion,  eq.  (118),  into  the  above  equation 
yields 


a<m)  (t)  =  J W(x)Zj  (x)  X Ik  (*Kk  -  E 

k=1 


X!k(x)vak 


k=l 


>dx 


=  J  W(x)Zi  (x)Xlk  (x)vakdx  -J  W(x)Zi  (x)E 

k=l 


ak 


k=l 


dx  .(121) 


Notice  that  the  second  term  is  the  product  of  piston  and  a  non-piston  Zernike  polynomial. 
Using  the  orthogonality  property  of  the  Zernike  polynomials  shows  that  the  second  term 
is  identically  zero.  Thus, 


74 


a|m)(t)  =  Jw(x)Zi(x)Xlk(x)vakdx 

k=l 


=  X[Jw(x)Zi(x)Il(x)dx]v,1 


T 

=  mi  Va 


(122) 


where 

=[Jw(x)Zj(x)I1(x)dx  J W(x)Z; (x)I2 (x)dx  •••  Jw(x)Zj(x)InA  (x)dxj  (123) 


and 

va=[val  va2  •••  vanJT.  (124) 

Therefore, 

a(2m)(t)  =  m>a 

ar(t)  =  mjva  ^a<m)(t)  =  MVa  (125) 

a(nT!(0  =  mn+lVa 

where  mjr  has  size  1  x  nA,  va  has  size  nAxl,  and  M  is  n  x  nA  where  n  is  the  number  of 
Zernike  modes  used  in  the  expansion  and  nA  is  the  number  of  actuators.  In  physical 
terms,  the  matrix  M  can  be  considered  the  steady-state  influence  function  matrix  since  it 
relates  control  voltages  to  the  steady-state  influence  of  the  mirror.  Complete  specification 
of  matrix  M  can  be  found  in  Appendix  B. 


5.3  Actuator  Dynamics 

Now  that  the  steady-state  behavior  of  the  mirror  has  been  considered,  it  is 
necessary  to  model  its  transient  behavior.  Piezoelectric  actuators,  used  in  the 
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construction  of  deformable  mirrors,  can  be  modeled  as  first-order,  linear,  time-invariant, 
deterministic  systems  [33].  Thus,  the  dynamics  equation  for  the  displacement  of  an 
actuator  can  be  modeled  by  the  scalar  equation 

8(t)  =  —  6(t)  +  bu(t)  (126) 

where  u(t)  is  the  command  signal  (voltage)  applied  to  the  actuator,  6(t)  is  the  path  length 
decrease,  and  Ta,  b  are  free  parameters  to  be  determined  by  empirical  data.  In  fact,  xa  is 
the  time  constant  of  the  actuators.  Under  the  assumption  that  all  actuators  behave  in  the 
same  manner,  nA  equations  identical  to  the  above  would  be  written.  As  discussed  earlier, 
the  states  in  this  case,  8(t),  are  the  2-way  decrease  in  path  length  since  the  light  is 
reflected  from  the  mirror. 

Now,  it  is  desired  to  formulate  the  dynamics  of  the  mirror  in  the  form 

•  (m) 

a  (t)  =  Fma(m)(t)  +  Bmu(t)  (127) 

where  a<m)  (t)  =  [a(2m)  (t)  a(,m)  (t)  •  •  •  a^j  (t)J  T  and  Fm,  Bm  are  the  dynamics  and 

control  input  matrices,  respectively.  Since  all  actuators  have  the  same  dynamics 
equation,  the  entire  mirror  exhibits  dynamic  behavior  with  a  time  constant  of  Ta.  Thus, 
the  mirror’s  contribution  to  the  phase  distortion  of  the  reflected  image  [33],  denoted  by 
<b(m)(x,t),  can  be  written  as 

.(m)  _] 

O  (x,  t)  =  — 0(m)(x,  t)  +b^(x)u(t)  (128) 

where  (x)  maps  control  voltages  to  the  rate  of  change  of  phase  at  spatial  location  x . 
Zernike  polynomials  can  also  be  used  to  write 
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0(m>(x,t)  =  Xa|m)(t)Zi(x). 


(129) 


Taking  the  time  derivative  of  the  above  equation  yields 


•  (m)  .(m) 

O  (x,t)  =  2>i  (t)Zi(x). 


(130) 


Substituting  eqs.  (129)  and  (130)  into  eq.  (128)  produces 

•  (m)  .(m)  •  (m) 

a2  (t)Z2(x)  +  a3  (t)Z3(x)+  •••  +  a„+i(t)Zn+1(x) 

=  —  [a(2m)(t)Z2(x)  +  a;;,)(t)Z1(x)+---+a^(t)Zn+1(x)]+bI(x)u(t).  (131) 

In  more  compact  notation,  eq.  (131)  becomes 

Xai  (t)Z,(x)  =  —  Xa|m)(t)Zi(x)  +  Xbj(x)uj(t)  (132) 

i=2  i=2  j=l 

where  bJ(x),uj(t)  are  the  jth  components  of  b^(x),u(t),  respectively.  Multiplying  the 
above  equation  by  W(x)Z2  (x)  and  integrating  over  the  area  of  the  mirror  yields 

n  +  l  .  (m)  _i  n+1  nA 

J  £ ai  (t)Z,  (x)W(x)Z2(x)dx  =— J  £ a,(m)  (t)Z,  (x)W(x)Z2  (x)dx  +  j  £  b f(x)u j(t)W(x)Z2  (x)dx  • 

i=2  Ta  i=2  j=l 

(133) 

Utilizing  the  orthogonality  property  of  Zemike  polynomials  produces 

.(m)  _] 

a2  (t)  = — a(m)(t)  +  r2Tu(t).  (134) 

" 

Performing  the  same  procedure  for  Z3(x),---,Zn+1(x)  gives 
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.(m)  _i 

a2  (t)  =  — a<m)(t)  +  r2Tu(t) 


•  (m) 


Xa 

-1 


a 3  (t)  =  -a"’(t)  +  r  u(t) 
x„  ' 


.(m)  _i 

an+i  (t)  =  — -a*™}  (t)  +  rj+,u(t) 

Tn 


(135) 


where  it  can  be  seen,  from  eqs.  (133)  and  (134),  that 

=  [| b, (x)W(x)Zj (x)dx  jb2(x)W(x)Z,(x)dx  •••  Jb„A (x)W(x)Zs (x)dx]  (136) 


and 


u(t)  =  [u,(t)  u2(t)  •••  un  (t)]T. 


(137) 


Thus, 


•  (m) 

a  (t)  =  DIAG 


vTay 


a(m)(t)  +  rTu(t) 


(138) 


where  DIAG(-l/xa)  is  a  diagonal  matrix  with  (-l/xa)  along  the  main  diagonal  and  zeros 


•  (m) 


everywhere  else.  At  steady-state,  a  (t)  =  0.  Therefore, 


0  =  DIAG 


^-0 


VXa/ 


aLm)(t)  +  r‘uss(t) 


(139) 


where  the  subscript  ‘ss’  implies  steady-state.  Rearranging  produces 

aLm)(t)  =  -DIAG(-xa)rTuss(t).  (140) 

Using  the  steady-state  information  developed  earlier,  that  is,  eq.  (125),  gives 

Muss(t)  =  -DIAG(-xJrTuss(t)  =»  rT  =  -DIAG^  jM .  (141) 


Hence,  the  mirror  model  becomes 
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•  (m) 

f-i) 

f-l] 

a  (t)  =  DIAG 

a(m)  (t)  + 

-DIAG 

Mu(t) 

U  7 

UJ 

=  Fma<m)(t)  +  Bmu(t) . 


(142) 


5.4  Complete  System  Model 

Thus  far,  two  separate  models  have  been  developed.  The  first  is  the  atmospheric 
model,  viz.,  the  plant,  which  is  a  linear,  stochastic,  dynamical  system  driven  by  Brownian 
motion.  This  model,  described  in  eq.  (116),  is  given  by 

da(a)(t)  =  Fa  a(a)(t)dt  +  Ga  dBa  (t)  (143) 

or,  in  less-rigorous,  white-noise  notation, 

aa>(t)  =  Faa(a)(t)  +  Gawa(t).  (144) 

Here,  the  state  vector  a(a)(t)  contains  the  coefficients  of  the  Zernike  polynomials  which 
expand  the  distortion  in  the  phase  of  the  inbound  (beacon)  and  outbound  (HEL) 
wavefronts.  The  second  is  the  deformable  mirror  model,  eq.  (142)  and  is  given  by 

•  (m) 

a  (t)  =  Fma<m)(t)  +  Bmu(t)  (145) 

where  a(m)(t)  contains  the  coefficients  of  the  Zernike  polynomials  which  expand  the  DM’s 
shape.  The  complete  system  model  is  formed  by  appending  these  two  systems: 


a(t)  = 


.(a) 

a  (t) 

•  (m) 

a  (t) 


Fa  0 
0  Fm 


a(a)(t) 

a(m)(t) 


0 

B„ 


u(t)  + 


Ga 

0 


wa(t) 


=  F a(t)  +  B  u(t)  +  Gwa(t) . 


(146) 
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An  output  equation  must  also  be  specified.  This  development  is  described  in  the  next 


section. 


5.5  Output  Equations 

Recall  that  the  state  vector  is  given  by 


a(t)  = 


1 

ft’ 

_ 1 

~a(a)(t)" 

__ 

a(0)(t) 

a<m)(t)_ 

1 

+-* 

? 

_ 1 

(147) 


The  quantities  that  can  be  measured  in  this  case  are  the  Zernike  coefficients  associated 
with  the  inbound  and  reflected  wavefronts.  Measurements  of  the  reflected  wavefront  are 
considered  only  if  the  metrological  laser  shown  in  Figures  2  and  6  is  used.  Considering 
the  inbound  wavefront  first,  a  matrix  selector  can  be  used  to  extract  the  relevant  Zernike 
coefficients,  that  is, 

a(i)(t)  =  [i  0  0]a(t)  =  E,a(t)  (148) 

where  I  and  0  are  n  by  n  identity  and  zero  matrices,  respectively,  and  E*  is  a  selection 
matrix.  It  is  desired  to  derive  the  measurement  equation  in  the  form 

zi(tj)  =  Hia(tj)  + Vj(tj)  (149) 


where  Vj(tj)  is  a  white,  Gaussian,  zero-mean,  discrete-time  process  with  variance 

E{vi(tj)vT(tk)}  =  Ri(tj)5jk.  (150) 

It  is  assumed  that  the  measurement  device  is  a  Hartmann  wavefront  sensor 
(HWFS)  [11,  23].  Basic  operation  of  this  sensor  is  as  follows:  the  HWFS  contains 
square  subapertures,  each  of  which  focuses  its  share  of  incident  light  onto  a  reticon 
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detector.  The  location  of  the  focused  light  on  the  detector  provides  a  measurement  of  the 
average  x-  and  y-tilts.  Thus,  the  HWFS  is  essentially  a  slope  sensor.  Hence,  phase 
distortion  in  the  inbound  wavefront’s  image  manifests  itself  in  a  set  of  subaperture  tilt 
measurements.  For  a  WFS  with  p  subapertures,  there  are  2p  outputs  (p  x-tilts  and  p  y- 
tilts).  It  is  further  assumed  that  the  HWFS  directly  outputs  these  slope  measurements. 

To  begin  the  derivation,  consider  the  phase  distortion  in  the  inbound  image  at 
discretized  times,  <F(1)(x,  t ; ),  where  x=  (X,Y)  defines  a  rectangular  position  on  the 


aperture  with  respect  to  the  aperture  center.  Since  the  HWFS  is  a  slope  sensor,  the  x-tilt 
prior  to  reaching  the  square  subaperture,  at  position  Xi,  Yi,  is 


x  -  tilt 


30(,)(x,tj) 


dY 


X,.Y, 


(151) 


Let  Xs,  Ys  define  the  center  of  the  sth  subaperture  with  respect  to  the  entire  sensor  array. 
Then,  the  average  x-tilt  going  into  the  sth  subaperture  is  [33] 

1  3®(X,Y,tj) 


_L  fs+~f 

A  Jv  J 


2 

A  •’y,-— Jx,-—  dY 
‘  2  2 


-dXdY 


(152) 


where  A  is  the  area  of  a  subaperture.  From  Peterson  and  Cho  [22],  the  output  of  a  HWFS 
is  a  similar  integral,  but  with  the  slopes  spatially  weighted.  In  particular,  for  the  x-tilt 
channel,  the  output  is 


s  7 

VA 
2 

where  L  is  a  constant  which  depends  on  wavelength,  lenslet  focal  length,  and  output 
scaling  and  vs  (tj)  is  an  element  from  the  vector  vs(tj).  Peterson  and  Cho  [22]  have  also 


rx*+ 


Va 


x - 


2 

VA 

2 


Wx  (X  -  Xs,  Y  -  Ys) 


3<D(X,Y,tj) 

dY 


dXdY  +  V:  (tT  (153) 

xs  ^ 
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derived  a  normalized  version  of  the  spatial  weighting  function,  which  for  the  x-tilt 
channel  can  be  expressed  as 


Wx(X  ,  Y  )  = 


2  In  2 


f 


21n2- 


1- 


V 


2Y 

Va" 


In 


V  Va 


A 


/ 


2Y 


A 


1+  /- 

v  Va  j 


In 


1  + 


2Y 


.  (154) 


The  inbound  wavefront’s  distorted  phase  image  can  be  written  as  a  linear 
combination  of  Zemike  polynomials,  that  is. 


0(i)(X,tj)  =  Xanil(tj)Zn+1(x). 

n 


(155) 


Taking  the  partial  derivative  of  0(1)(x,tj)  with  respect  to  Y,  as  described  in  eq.  (151), 
produces 


ao(i)(x,tj) 

3Y 


=  Sa-Vtj) 


azn+,(x) 

3Y 


(156) 


Using  eq.  (156)  in  eq.  (153)  yields 


y - X. -  n 

s  2  2  n 


3Zn+,(x) 
3  Y 


l-x 

A^ 


fY‘+—  fA‘+—  w  v  Y  Y^Zn+,(x) 

L  VaJx  Va  Wx(X-Xs,Y-Ys) 

s  2  s  2 


dXdYa^tp  +  v^tj) 


L 

A 


nT  o(i) 

Px  a 


(tj)  +  vix  (tj) 


Eia(tj)  +  vi  (tj) 

J  xs  J 


(157) 


where 
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(158) 


[ffw„3Z’®dXdY 

rr  az3(x) 

ffw  ^Z"+l  dXdyl 

JJ  *  ay 

JJ  *  ay 

JJ  ay 

and  Wx  =  Wx(X-Xs,  Y-Ys).  In  order  to  specify  measurements  of  the  inbound  wavefront 
completely,  expressions  similar  to  those  in  eq.  (157)  would  have  to  be  written  for  y-tilts 
and  for  each  of  the  p  subapertures.  Putting  these  equations  into  matrix  form  gives 

zi(tj)  =  -t-NEia(tj)  +  vi(tj) 

=  H,a(tj)  +  V|(tj)  (159) 

where  H;  =— NEj  and  N  is  a  2p  x  n  matrix  with  the  following  form: 

A 

r  -lT 

N=  px,  Px2  •••pxp  Py|  Py2  •••Pyp  •  (160) 

More  information  on  matrix  N  can  be  found  in  Appendix  B. 

Proceeding  in  a  similar  manner  for  measurements  of  the  reflected  wavefront,  the 
selection  matrix  is 

a(m,(t)  =  [0  0  l]a(t)=Ema(t)  (161) 

where  I  and  0  are  n  by  n  identity  and  zero  matrices.  The  measurement  equation  in  this 
case  becomes 

zm(tj)  =  Hma(tj)  +  vm(tj)  (162) 

where  vm(tj)  is  a  white,  Gaussian,  zero-mean,  discrete-time  process  with  variance 

E{vm(tj)v^,(tk)}  =  Rni(tj)5jk  (163) 

and  Vi(tj)  and  vm(tj)  are  independent  random  vectors.  In  this  case,  the  x-tilt  output  of  the 
HWFS  is 
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-.<v=LxfMw.(x-x--Y-)2a»,p^— 

^  *s  2  2  n 


3Zn+l(x) 
dY 


=  yPl  Ema(tj)  +  vmSs(tj). 


Writing  equations  for  each  subaperture  and  for  the  y-tilts  produces 


Zm(tj)  =  ANEn,a(tj)  +  Vn,(tj) 


Hma(tj)  +  vm(tJ) 


(164) 


(165) 


where  Hm  = — NEm. 

m  A 


In  summary,  the  following  dynamics  model  and  output  equation  have  been 
developed: 


a(t)  = 


.00 

a  (t) 

.  (ni) 

a  (t) 


Fa  0 
0  F„, 


a<a)  (t) 


a(m)  (t) 


Br 


u(t)  - 


G. 

0 


w  (t)  =  F  a(t)+B  u(t)  +  Gw,(t)  >  (166) 


Zj(tj)  =  H,a(tj)  +  Vj(tj) , 
zni(lj) =  Hma(tj)  +  vm(tj) , 


Z(tj)  = 


’zi(tj)’ 

"Hi  " 

a(tj)  + 

"vi(tj)" 

_Zm(t  j>_ 

_Hm 

_Vm  (t  j )_ 

=  Ha(tj)  +  v(tj) , 


(167) 

(168) 

(169) 


where  wa(t)  is  a  zero-mean,  white,  Gaussian  noise  with  E{wa  (t) wa  (t  +  x)}  =  1 8(x) ,  Vi(tj) 
is  a  zero-mean,  white,  Gaussian,  discrete -time  process  with  E|vj(tj)vf  (tk)|  =  Rj(t  |)8jk , 
vm(tj)  is  a  zero-mean,  white,  Gaussian,  discrete-time  process  with 
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E{vm(tj)v:(tk)}  =  Rro(tj)8jk,  H  = 


Hi 

HL 


,  and  v(tj)  = 


VjCtj) 

Vm(tj) 


.  Also,  v(tj)  is  a  zero- 


mean,  white,  Gaussian,  discrete-time  process  with  variance 

YRi(tj)  o  > 
o  Rm(tj) 


E{v(tj)vT(tk)}  = 


E{vj(t  j)vs  (tk )}  E{v.(tj)vm(tk)} 
|_E{vm(tj)Vj(tk)}  E{vm(tj)vm(tk)}J 


=  R(tj)  j  =  k 


"0  0 
[o  0, 


(170) 


Equation  (166)  represents  the  underlying  dynamic  system  for  the  Kalman  filter. 


5.6  Discrete-Time  System 

The  system  in  eq.  (166)  is  a  continuous-time  model  of  the  atmosphere  and 
deformable  mirror.  Since  the  Kalman  filter  and  controller  to  be  developed  in  the  next 
chapter  will  be  implemented  in  a  digital  computer,  it  is  necessary  to  convert  eq.  (166)  to 
an  equivalent  discrete-time  system  [16].  To  perform  this  conversion,  consider  the 
solution  to  eq.  (166)  which  can  be  written  as 


a(tJ„)  =  ®(tw,tJ)a(tJ)  + 

1  j+1 

J<I>(tj+I,x)Bu(x)dx 

_  *j 

+ 

lj+l 

J<D(tj+,,T)Gwa(x)dx 

=  4>(tj«.tj)a(tJ)  + 

lj+l 

}o(tj+1,x)Bdt 

U(tj)  + 

lJ+l 

J  <I>(tj+1,x)Gwa(x)dx 

=  4>(tj+1,tj)a(tj)  +  Bd(tJ)u(tj)  +  wd(tj).  (171) 

Note  that  the  second  expression  in  eq.  (171)  can  be  written  since  it  is  assumed  that  the 
control,  u(t),  is  held  constant  over  each  sample  period.  Also,  in  eq.  (171),  <&(tj+i,tj)  is  the 
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state  transition  matrix,  Bd(tj)  =  J <b(tj+l  ,x)Bdx  ,  and  wd(tj)  is  a  white,  Gaussian, 

t  j 

discrete-time  stochastic  process  with  statistics  duplicating  those  of 

1  i+J 

|  <b(tj+l  ,x)Gwa(x)dx.  In  other  words,  wd(tj)  has  statistics 

‘j 

E{wd(ti)}  =  0 

lj+l 

E{wd(tj)Wj(tj)}  =  Qd(tj)=  J<I>(tj+1,x)GIGT<l>T(tj+1,x)dx  (172) 

•i 

e{w<i (t j)wd (tk )}  =  0  j*k. 

Since  the  complete  system  model  is  linear  and  time-invariant,  <J>(tj+i,tj)  reduces  to  a 
function  of  the  time  difference  At  =  tj+i  -  tj.  Therefore, 

^(t  j+,  ,t  j)  =  0(tj+l  -  t  j)  =  eF<,jtl_,j) .  (173) 

Appendix  B  contains  a  complete  explanation  of  the  determination  of  the  discrete-time 
system  matrices  for  simulation  purposes. 

5.7  Summary 

At  this  point,  all  of  the  required  modeling  has  been  performed.  In  particular, 
atmospheric  and  mirror  models  have  been  formed,  output  equations  have  been  specified, 
and  an  equivalent  discrete-time  system  has  been  developed.  Now,  these  models  can  be 
embedded  in  a  Kalman  filter  so  that  estimation  of  the  outbound  wavefront  Zernike 
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expansion  coefficients  can  be  performed.  Using  these  estimates,  a  linear  quadratic 
Gaussian  controller  will  apply  the  estimated  conjugate  phase  to  the  mirror.  The  next 
chapter  addresses  design  of  the  filter  and  controller. 
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VI.  Kalman  Filter  and  LQG  Controller 


6.1  Introduction 

The  objective  of  this  research  is  to  design  a  control  system  which  can  provide 
commands  to  the  deformable  mirror  to  cancel  the  effects  of  atmospheric  turbulence  in  the 
highly  dynamic  scenario  of  the  ABL.  The  atmospheric  model,  mirror  model,  and 
measurement  equations  have  been  discussed  and  are  given  by  eqs.  (166)  -  (170). 


6.2  Kalman  Filter 


A  Kalman  filter  is  used  to  estimate  the  state  vector  a(tj).  Recall  that  a(tj)  is 
composed  of  the  inbound,  outbound,  and  mirror  Zernike  polynomial  expansion 


coefficients,  i.e.,  a(t)  = 


{a">(t)}T  {a'”'(0}T  {a'-O)}1 


IT 


.  The  Kalman  filter  accepts 


measurements,  z(tj),  and  outputs  estimates  of  the  system  state,  a(t ;)  ,  where  the 

superscript  ‘A’  is  used  to  denote  an  estimate.  The  filter  performs  these  estimates  using  an 
internal  model,  eq.  (171),  propagating  the  state  estimate  and  state  estimation  error 
covariance,  and  then  appropriately  incorporating  the  measurements,  eq.  (167).  The 
filter’s  state  estimation  error  covariance,  denoted  P|{t),  is  defined  by 

P,  (t)  =  E{[a(t)  -  a(t)]  [a(t)-  a(t)]T|Z(t)}  (174) 


where  the  expectation  is  conditioned  on  the  measurement  history  Z(t). 

Operation  of  the  Kalman  filter  can  be  divided  into  two  processes:  propagation  and 
measurement  update.  Propagation  is  the  change  in  state  estimate  and  filter  covariance  as 
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time  is  moved  forward  from  tj.i  to  tj.  Measurement  update  is  the  change  in  state  estimate 
and  filter  covariance  as  measurements  are  processed.  The  governing  equations  for  these 
two  cycles  are  [16]: 


a(t: )  =  <D(t .  - 1 H  )a(tj_, )  +  Bd (t H )u(tH ) , 

(175) 

Pf(t:)  =  0(tj-tH)Pf(t;i)®T(tj-tH)  +  Qd(tj), 

(176) 

K(tj)  =  Pf  (t:)HT[HPf  (t:)HT  +  R(t .)[' , 

(177) 

a(tp  =  a(tT)  +  K(tj)[zj-Ha(tT)], 

(178) 

Pf(t|)  =  Pf(t:)-K(tj)HPf(t:) 

(179) 

where 

a(tj.,)  = 

state  estimate  at  the  start  of  the  tj.i  to  tj  propagation  cycle  and  after  the 
measurement  update  at  tj.i. 

Pf(tn)  = 

filter  covariance  at  the  start  of  the  tj.t  to  tj  propagation  cycle  and  after  the 
measurement  update  at  tj.j, 

a(t-)  = 

state  estimate  at  the  end  of  the  tj_i  to  tj  propagation  cycle  and  before  the 
measurement  update  at  tj, 

Pf(tJ)  = 

filter  covariance  at  the  end  of  the  tj_i  to  tj  propagation  cycle  and  before  the 
measurement  update  at  tj, 

K(tj)  = 

Kalman  filter  gain, 

a(tp  = 

state  estimate  at  the  end  of  the  tj.i  to  tj  propagation  cycle  and  after  the 
measurement  update  at  tj, 

Pf(tp  = 

filter  covariance  at  the  end  of  the  tj.i  to  tj  propagation  cycle  and  after  the 
measurement  update  at  tj, 
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and 


Zj  =  measurement  realization  at  the  update  time. 

Initial  conditions  on  the  filter,  that  is,  a(t0)  and  Pf(to),  must  also  be  specified.  Therefore, 
the  initial  conditions  and  eqs.  (166)  -  (171),  (175)  -  (179)  fully  describe  the  Kalman  filter. 

6.3  LQG  Controller 

In  this  section,  an  LQG  (Linear  dynamic  system  model,  Quadratic  cost,  and 
Gaussian  noise  inputs)  controller  [17]  will  be  developed.  The  general  idea  here  is  that  the 
behavior  of  some  underlying  dynamic  system  is  to  be  affected  by  control  inputs  in  such  a 
way  that  the  controlled  variables  exhibit  some  desired  characteristics  [17].  The  controller 
will  in  fact  be  a  regulator,  which  is  based  on  the  assumption  that  there  are  costs 
associated  with  nonzero  controlled  variables  and  nonzero  controls.  The  cost  function  is 
quadratic  in  nature,  that  is,  the  cost  is  proportional  to  the  weighted  squares  of  the 
controlled  variables  and  inputs. 

For  the  ABL,  the  inbound  wavefront’s  phase  is  measured  to  determine  the 
turbulence-induced  wavefront  aberrations.  It  is  then  desired  to  apply  the  estimated 
conjugate  (opposite)  phase  to  the  mirror  so  that  the  HEL’s  wavefront  phase  perturbation 
will  be  minimal  upon  reaching  the  target,  viz.,  the  Strehl  ratio  at  the  target  is  maximized. 
Therefore,  the  controlled  variables  are 

y(tj)  =  [0  -I  l]a(tj)  =  Eca(tj)  =  -a(0)(tj)  +  a(m)(tj)  (180) 

where  the  size  of  y(tj)  is  n  by  1 . 
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A  useful  form  of  quadratic  cost  function  is  [17] 


+  uT(t,)U(tj)u(tJ)]  +  V(t 


LN+1 


)Yfy(tN+1) 


(181) 


where  Y(tj)  V  tj  and  Yf  are  selected  as  n  by  n  (n  =  2  for  tilts)  positive  definite  matrices  and 
U(tj)  is  a  positive  definite  matrix  of  dimension  equal  to  the  number  of  actuators.  In  the 
simulations  discussed  in  Chapter  7,  the  quadratic  cost  function  matrices  were  selected  as 


Y(tj)  =  Yf 


5  0 
0  1 


,U(tj)  =  05I. 


(182) 


These  values  were  selected  after  many  iterations  to  provide  satisfactory  state  regulation 
without  expending  too  much  control  energy.  Recall  that  the  controlled  variables  are 

y(t j)  =  [-a^0) (t)  +  aS,"0 (t)  -a^ (t)  +  (t)]T .  Since  dynamic  motion  occurs  only  in  the  x- 

direction,  it  is  expected  that  there  will  be  more  error  in  the  x-tilt  controlled  variables  as 
compared  to  the  y-tilt  variables.  Hence,  entry  (1,1)  of  Y(tj)  and  Yf  was  set  to  5  to  place 
more  emphasis  on  regulation  of  the  x-tilt  controlled  variables. 

Equation  (181)  does  not  include  any  cross-terms  between  the  controlled  variables 
and  controls.  Cross-terms  are  useful  when  it  is  desired  to  exert  control  influence  on  the 
state  vector  over  the  entire  sample  period,  not  just  at  the  sample  points.  However,  since 
the  sample  period  is  small  for  this  system,  it  is  expected  that  cross-terms  would  be  of 
negligible  benefit  and  thus  will  not  be  necessary.  Using  the  expression  for  the  controlled 
variables,  eq.  (180),  in  eq.  (181),  the  cost  function  becomes 

J  =  E|^^[aT(tj)A(tj)a(tj)  +  uT(tj)U(ti)u(tj)]  +  ^-aT(tN+1)Afa(tN+1)|  (183) 
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where  A(tj)  =  EcTY(tj)Ec  and  Af  =  EcTYf  Ec.  It  is  desired  to  determine  the  control 
function  u*(«)  which  minimizes  the  cost  function  given  in  eq.  (183).  It  can  be  shown  that 
the  cost  minimizing  controller  [17]  is 

u*(tj)  =  -G*(tj)a(tJ)  (184) 

where  Gc*(tj)  is  the  controller  gain.  This  gain  can  be  calculated  from  the  solution  to 

j)  =  [^(tj)  +  Bd(tj)Kc(tj+1  )Bd(tj)]"'[Bj(tj)Kr(tj+l)0(tj+1,tJ)]  (185) 

where  Kc(t,)  satisfies  the  backward  Riccati  equation  [17] 

K,(tj)  =  A(tJ)  +  4>T(tJ„,tJ)KI(tj„)®(tH,tJ)-[<l>T(tJtl,tj)Kc(tH)BJ(tj)]G;(ti) 

=  A(t|)  +  4>T(tH,tJ)K,(ej„)[<l.(tJ.l.tJ)-Ba(ti)G;(tJ)]  (186) 

solved  backwards  from  the  terminal  condition 

Kc(tN+1)  =  Af .  (187) 

The  filter  and  associated  measurements  devices  are  denoted  by  Kalman  filter  1 , 

Kalman  filter  2  and  WFS  1,  WFS  2  in  Figure  25.  The  input  to  Kalman  filter  1  is  zj(tj), 
that  is,  the  measured  Zernike  coefficients  of  the  inbound  wavefront,  while  the  input  to 
Kalman  filter  2  is  zm(tj),  that  is,  the  measured  Zemike  coefficients  of  the  reflected 
wavefront.  Measurements  zm(tj)  are  included  in  this  figure  in  case  the  metrological  laser 
is  used.  Currently,  a  metrological  laser  is  not  used  in  the  ABL  and  thus  these 
measurements  are  not  incorporated  into  the  simulations.  The  blocks  ACT  and  DM  in 
Figure  25  represent  the  actuator  dynamics  and  deformable  mirror,  respectively.  Also,  Gc* 
and  u*  represent  the  LQG  controller  gain  and  the  cost  minimizing  control  function.  The 
atmospheric  and  DM  models  contain  no  coupling  between  each  other,  therefore,  the 


92 


single  augmented  Kalman  filter  (atmospheric  and  DM  models)  can  be  decomposed  into 


the  two  decoupled  filters,  Kalman  filter  1  and  Kalman  filter  2,  displayed  in  Figure  25. 
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Figure  25.  Adaptive  optics  control  system. 


The  dimensionality  of  the  variables  associated  with  the  control  system  is  given  in  Tables 
2a  and  2b.  In  these  tables,  n  is  the  number  of  Zernike  modes  used  in  the  wavefront 
expansions,  p  is  the  number  of  subapertures  in  the  Hartmann  wavefront  sensor,  and  nA  is 
the  number  of  actuators. 
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Matrix/Vector 

Dimension 

a(t) 

3n-by- 1 

aw(t) 

2n-by- 1 

a(m)(t) 

n-by- 1 

a(i)(t) 

n-by- 1 

a(0)(t) 

n-by- 1 

m 

nA-by- 1 

M 

n-by-nA 

va 

nA-by- 1 

Tj 

nA-by-l 

r 

nA-by-n 

DIAG(1  /x) 

n-by-n 

Ej  and  Em 

n-by-3n 

Fa 

2n-by-2n 

Fm 

n-by-n 

F 

3n-by-3n 

<E>(ti+i,ti) 

3n-by-3n 

Bm 

n-by-nA 

B 

3n-by-nA 

u(t) 

nA-by- 1 

Ga 

2n-by-2 

G 

3n-by-2 

Wa(t) 

2-by-l 

Zi(tj) 

2p-by- 1 

Hi 

2p-by-3n 

Ri(ti) 

2p-by-2p 

Zm(tj) 

2p-by- 1 

Vi(tj) 

2p-by-l 

Table  2a.  Variable  dimensions. 
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Matrix/V  ector 

Dimension 

Hm 

2p-by-3n 

Vm(tj) 

2p-by-l 

Rm(tj) 

2p-by-2p 

Z(tj) 

4p-by-l 

H 

4p-by-3n 

V(tj) 

4p-by- 1 

R(ti) 

4p-by-4p 

N 

2p-by-n 

PxT,  PvT 

1-by-n 

Bd(ti) 

3n-by-nA 

Wd(tj) 

3n-by-l 

Qd(ti) 

3n-by-3n 

a(t;_,),a(t:),a(tp 

3n-by-l 

Pf(tn),  Pf(t:),  Pf(t|) 

3n-by-3n 

K(tj) 

3n-by-4p 

Zi 

4p-by-l 

Ec 

n-by-3n 

y(ti) 

n-by-1 

U(tj) 

nA-by-nA 

A(tj),  Af 

3n-by-3n 

U  (tj) 

nA-by-l 

Gc*(ti) 

nA-by-3n 

U(ti) 

nA-by-nA 

Y(tj),  Yf 

n-by-n 

Y(tN+i) 

n-by-n 

Ke(ti) 

3n-by-3n 

Kc(tN+l) 

3n-by-3n 

Table  2b.  Variable  dimensions. 
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6.4  Summary 

A  Kalman  filter  and  LQG  controller  were  designed  in  the  last  few  sections.  The 
filter  produces  estimates  of  the  state  variables  while  the  controller  is  used  to  drive  the 
controlled  variables  to  zero.  That  is,  the  controller  generates  gains  which  are  used  to 
provide  commands  to  the  DM. 

At  this  point,  all  of  the  analysis  has  been  set  into  place.  Therefore,  the  complete 
system,  that  is,  atmospheric  and  mirror  models,  output  equations,  filter  and  controller  can 
be  simulated  to  evaluate  performance.  To  accomplish  this  simulation,  a  phase  screen, 
which  will  be  discussed  in  Chapter  7,  will  be  generated  to  provide  the  measured  quantity. 
Simulation  results  involve  the  Root  Mean  Square  (RMS)  phase  distortions,  Strehl  ratios, 
phase  screen  and  estimated  quantities,  and  Monte  Carlo  statistics.  The  next  chapter 
addresses  the  simulations  and  analyzes  the  results. 
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VII.  Simulations 


7.1  Introduction 

Now  that  all  of  the  modeling  has  been  performed,  the  entire  system  can  be 
simulated  to  evaluate  its  performance.  Before  the  simulations  can  be  run,  a  few  other 
items  must  be  covered,  including  performance  metrics,  the  phase  screen,  and 
development  of  the  simulated  correlation  functions. 

7.2  Performance  Metrics 

In  this  section,  the  metrics  used  to  evaluate  system  performance  are  delineated. 
These  quantities  provide  a  measure  of  the  performance  of  the  controller.  There  are  two 
metrics  that  will  be  used  in  this  work.  Imaging  performance  is  specified  by  a  quantity 
termed  the  Strehl  Ratio  (SR)  and  the  Root  Mean  Square  (RMS)  phase  distortion.  The 
Strehl  ratio  is  defined  as  the  irradiance  with  wavefront  aberrations  present  divided  by  the 
diffraction-limited  irradiance  on  the  optical  axis  [1],  A  useful  estimate  of  the  Strehl  ratio 
can  be  made  if  the  aperture  averaged,  squared,  residual  wavefront  phase,  <f>^es(t) ,  is  less 
than  (27E  /  10)2  rad2  [23].  In  this  case,  the  Strehl  ratio  becomes 

SR  ~  exp{-C>2es(t)} .  (188) 

From  eq.  (188),  it  is  obvious  that  if  02es(t)  =  0,  that  is,  perfect  correction  of  the 

wavefront  phase  deformations  is  performed,  then  SR  =  1.  This  is  the  upper  bound  on  the 
Strehl  ratio.  In  most  applications,  a  Strehl  ratio  of  0.6  or  higher  is  adequate.  However,  it 
should  be  noted  that  the  Strehl  ratio  is  dependent  on  the  seeing  conditions  (clear, 
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overcast,  precipitation,  etc.),  AO  hardware,  and  other  factors.  Hence,  a  measure  of  a 
“good”  Strehl  ratio  is  difficult  to  present. 

There  are  three  Strehl  ratios  that  are  of  interest  in  this  work.  The  first  is  the  Strehl 
ratio  of  the  outbound  wavefront  which  provides  a  measure  of  the  uncompensated  system 
performance.  The  second  is  the  Strehl  ratio  of  the  corrected  wavefront  using  simple 
phase  conjugation.  Simple  phase  conjugation  implies  that  measurements  of  the  inbound 
wavefront’s  phase  deformations  (x-  and  y-tilts)  are  directly,  negatively,  applied  to  the 
DM.  In  this  case,  there  are  no  correlation  functions  to  determine.  It  is  still  a  closed-loop 
system;  however,  anisoplanatism  is  not  taken  into  account.  This  is  the  most  rudimentary 
phase  conjugation  scheme  available  and  is  currently  used  in  most  AO  applications. 

Simple  phase  conjugation  provides  an  algorithm  on  which  to  base  the  performance  of  the 
AO  tilt  compensation  algorithm  developed  in  this  work.  The  last  Strehl  ratio  is  that  of  the 
compensated  wavefront  obtained  from  the  work  in  previous  chapters.  In  this  case,  the 
necessary  correlation  functions  are  calculated,  estimates  of  the  outbound  wavefront 
Zernike  phase  expansion  coefficients  are  obtained,  and  the  LQG  controller  attempts  to 
drive  the  controlled  variables  to  zero.  Comparing  these  three  Strehl  ratios,  which  involve 
increasing  system  complexity,  provides  a  measure  of  the  performance  enhancement  that 
can  be  obtained  using  tilt  compensation. 

Since  the  wavefronts  have  been  expanded  using  Zernike  polynomials,  the  Strehl 
ratios  take  on  a  form  involving  the  Zernike  expansion  coefficients.  In  particular,  for  the 
fully  compensated  wavefront,  the  residual  wavefront  phase  becomes 

<I)reS(x,t)  =  -0(0)(x,t)  +  Cl>(nl>(x,t)  (189) 
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where  <I>(o)(x,t),<I>(rn)(x,t)  are  the  outbound  wavefront  phase  and  the  phase  correction 
applied  by  the  DM,  respectively.  The  aperture  averaged,  squared,  residual  wavefront 
phase  is 

®L(l)  =  Jw(x)OL(>i,t)<S  (190) 

where  <lJ^es(x,t)  is  the  squared,  residual  wavefront  phase.  Using  the  Zemike  polynomial 
expansions,  that  is, 

•“Ot.O-Z^Mz.OO  -  =  (19» 

j  k 

in  eq. (190)  produces 

Up  to  this  point,  only  x-  and  y-tilts  have  been  considered  in  wavefront  expansions. 
However,  in  terms  of  computing  Strehl  ratios,  it  is  advantageous  to  use  more  Zernike 
modes  in  the  inbound,  and  therefore,  outbound,  wavefront  expansions.  If  the  inbound 
wavefront  was  described  entirely  by  tilts,  then  the  Strehl  ratios  would  be  unusually  high, 
since  tilts  are  removed  from  the  corrected  wavefront.  For  example,  assume  that  the 
inbound  wavefront  phase  is  only  corrupted  by  tilts.  Then,  after  compensating  and 
applying  the  conjugate  phase  to  the  mirror,  the  residual  wavefront  phase  would  be  near 
zero,  assuming  the  filter/controller  performed  well.  In  this  case,  the  Strehl  ratio  would  be 
near  one,  that  is,  the  diffraction  limit.  Unfortunately,  this  type  of  AO  performance  is  not 
feasible.  To  rectify  this  situation,  more  Zernike  modes  can  be  used  to  expand  the  inbound 
and  outbound  wavefronts.  In  practice,  an  accurate  representation  of  these  actual 


4>LW  =  jw(x) 


X a<« (tJZ, (5E)  +  S at"’  (t)Z„  (X)  dx  .  (192) 

J  k 
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wavefronts  would  require  an  infinite  number  of  Zemike  modes.  Obviously,  the 
expansion  must  be  truncated  at  some  point;  however,  tilts  are  not  the  only  components  of 
atmospherically  distorted  wavefront  phases.  In  fact,  five  modes,  x-  and  y-tilts,  focus,  and 
orthogonal  components  of  astigmatism,  were  used  to  provide  a  more  accurate 
representation  of  the  actual  inbound  and  outbound  wavefronts.  The  AO  scheme 
developed  in  this  work  only  compensates  for  the  tilt  modes,  therefore,  the  outbound 
wavefront,  after  compensation,  is  still  corrupted  by  the  remaining  three  modes.  In  this 
way,  more  accurate  Strehl  ratios  will  be  obtained. 

For  n  Zernike  modes  (n  =  5  for  tilts,  coma,  astigmatisms),  eq.  (192)  becomes 


4>LW  =  f W(x)[-a“(t)Z, (x)  —  - a<»  (t)Z„, (x)  +  a‘« (t)Z, (*)  +a™ (t)Z,(x)f  dx .  (193) 

Notice  that  tilt  compensation  results  in  including  only  a2(m)(t)  and  a3(m)(t).  Expanding  the 
term  in  brackets  in  eq.  (193)  and  using  the  orthonormality  property  of  the  Zernike 
polynomials  produces 


.  (194) 


Therefore,  the  third  Strehl  ratio  becomes 


SR3  - (I)  +  (t)]2  +  (t)  +  (t)]‘  +  (t)}2  +-  +  {a£>(t)}  jj}  -  (195) 

For  the  Strehl  ratio  of  the  outbound  wavefront,  the  residual  wavefront  phase  is 

®«.M  =  4>“(x,t)  =  2>S0)(t)Zi(x) .  (196) 


Then, 


too 


<s(t)  =  Jw(x) 


E^WZjOO 


dx  =  J  W(x)[a<0)  (t)Z2  (x)  +  •  •  •  +  a^j  (t)Zn+1  (x)f  dx .  ( 1 97) 


Applying  the  orthonormality  property  again  produces 


4>L(t)  =  {a?>(t)}2+{a“(t)}!+...+{a«„"»(t)}2. 


(198) 


Hence,  the  first  Strehl  ratio  becomes 


SR,  ~  exp<(-l) 


{a'2">(,)}2+{a<->(.)}2+-+{a«(t)}! 


(199) 


For  the  simple  phase  conjugation  Strehl  ratio,  a  little  more  analysis  is  needed.  In 
the  simple  conjugation  case,  what  is  measured,  that  is,  the  set  of  inbound  wavefront 
Zernike  coefficients,  is  negatively  applied  to  the  DM.  Therefore,  it  is  desired  to  have 

a(m)(t)  =  -a(i)(t) .  (200) 

In  eq.  (200),  it  is  assumed  that  only  tilt  compensation  is  performed.  Thus,  even  though 
five  Zernike  modes  are  used  to  represent  the  inbound  wavefront,  only  the  negatives  of  the 
tilt  modes  are  applied  to  the  mirror.  For  the  simple  phase  conjugation  case,  the  Kalman 
filter  developed  in  Chapter  6  is  not  needed.  Recall  that  the  filter  was  used  to  estimate  the 
outbound  wavefront  Zernike  coefficients.  With  simple  phase  conjugation,  no  outbound 
wavefront  estimates  are  required  since  only  the  set  of  inbound  wavefront  Zernike 
coefficients  is  applied  to  the  mirror.  Therefore,  a  weighted  least  squares  algorithm  can  be 
used  to  estimate  the  inbound  wavefront  Zernike  coefficients  from  the  noisy  measurements 
as  described  in  that  which  follows.  Of  course,  a  Kalman  filter  could  also  be  built  to 
estimate  the  inbound  wavefront  Zernike  coefficients,  but  the  performance  benefit  over  the 
least  squares  estimate  probably  would  not  be  substantial. 
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Unfortunately,  the  inbound  coefficients  are  not  directly  measured.  Instead,  recall 
that  the  measured  quantity  is 

zi(tJ)  =  H,a,i>(tj)  +  v,(tj).  (201) 

Therefore,  it  is  necessary  to  generate  an  estimate,  a0)(tj) ,  of  a(,)(tj) .  Hence  the  value 
of  a(0(tj)  which  minimizes  the  weighted  sum  of  squares  of  the  components  of  the  vector 

[z,(.J)-H,a">(tj)]  (202) 


is  desired.  Stated  in  another  fashion,  the  a0>  (t  j ) ,  which  minimizes  the  cost  function 


H;a 


'‘’(•jlfw-M'j) 


(203) 


must  be  evaluated.  In  eq.  (203),  the  subscript  ‘sc’  implies  simple  conjugation  and  Wsc  is 
the  weighting  matrix.  Minimization  of  the  cost  function  in  eq.  (203)  is  accomplished  by 
satisfying  the  following  criteria: 


3J. 


5jsc 

V 

II 

Dl 

1 

3if(t 

a=a 

=  0T 


32  J 


(204) 


>0 


where  a^ftjj.a^tj)  are  estimates  of  the  x-  and  y-tilts  of  the  inbound  wavefront. 
Performing  the  differentiation  in  eq.  (204)  yields 

a">(tJ)  =  [H^W„Hi]‘lH^W„z1(tJ).  (205) 


Therefore,  using  eq.  (200),  the  desired  phase  conjugation  scheme  becomes 

a<m>  (t  j) =  -**(i)  (t  j)  =  -[hT  W^Hj  ]”'  HT  WscZj  (t  j ) .  (206) 
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The  Strehl  ratio  for  this  case  is  the  same  as  that  given  in  eq.  (195)  with  a(m)(t)  replaced  by 
the  expression  in  eq.  (206). 

The  second  performance  metric  that  will  be  used  is  the  RMS  phase  distortion. 

This  quantity  can  be  expressed  as 

4>RMs(t)  =  1/jw(x)<I>2(x,t)dx.  (207) 


For  the  outbound  wavefront, 


Cs(t  )  =  Jjw(x) 


-|2 


dx 


=  1/|w(x)[a“(x)Z!(x)+-+a“1(x)Z„.,(x)]2dx  .  (208) 

Expanding  the  term  in  brackets  and  once  again  using  the  orthonormality  property  of  the 
Zernike  polynomials  gives 

4>Ss(t)  =  1/{ar(t)}2+-+{al">(0}2  =1|SbrW}2  •  (209) 

For  the  corrected  image, 

0(co,rected)(-jt)  =  ^_a(jO)(t)Zj(-)  +  ^a(m)(t)Zk(-)  (210) 

j  k 

Substituting  this  expression  into  eq.  (207)  yields 

2>“(<)z,(x)  + 

j 

Performing  some  simplifications,  as  detailed  in  Appendix  B,  produces  the  desired  result, 
that  is. 


I4m,(t)zk(x) 


dx  . 


(211) 


^(corrected) 

M  RMS  \l)~ 


Jw(x) 
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(corrected) 

^RMS 


This  then  provides  a  set  of  rather  simple  equations  to  be  used  as  the  metrics  in  evaluating 
system  performance.  The  next  section  discusses  the  simulations  in  more  detail  along  with 
a  description  of  the  phase  screen. 

7.3  Simulation  Methodology 

A  digital  simulation  of  the  adaptive  optics  system  was  desired  in  order  to  evaluate 
the  effectiveness  of  the  LQG  controller.  This  simulation  was  performed  using  code 
written  in  MATLAB  [15].  To  begin  the  simulation,  the  time  history  of  the  covariance 
matrices  was  developed.  This  was  accomplished  in  the  following  manner:  given  the 
form  of  the  system,  input,  and  covariance  matrices,  F,  G,  and  P(t), 

Pn  (t)  0  pI3(t)  0  f,,  0  f13  0  g,  0 

P(t)=  °  P“(,)  °  P“(,)  ,F=  °  f“  °  f“  ,G=  °  82  ,(213) 

P 1 3  ( t )  o  P33(t)  0  f3l  0  f33  0  gj  0 

0  p24(t)  0  P44(t)_  _  0  f42  0  fw_  _  0  g4_ 

form  the  system  of  equations  P(t)  =  FP(t)  +  P(t)FT  +  GGT  .  Performing  the  algebra 
required  by  the  differential  equation  discussed  above  provides 

F(t)  =[COLl  COL2  COL 3  COL4]  (214) 

where 

COLl  =  [2fllpI1(t)  +  2f13p1J(t)  +  gf  0  f31pll(t)  +  (fll  +  f33)pl3(t)  +  f|3P33(t)  Of 
COL2  =  [o  2f22p22 (t)  +  2f24p24 (t)  +  g2  0  f42p22(t)  +  (f22 +f44)p24(t)  +  f24p44(t)]T  ^215) 
COL3  =  [f3lpn(t)  +  (fM  +  f33)p  l3  (t)  +  fJ3p33(t)  0  2f33p33(t)  +  2f31p31  (t)  -t-  g3  Of 
COL4  =  [o  f42p22(t)  +  (f22 +f44)P24(t)  +  f24P44(t)  0  2f44p44(t)  +2f42p42(t)  +  g4f . 
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Now,  form  the  linear  system 


Pn  P 13  P22  P24  P33  P44J  (0  —  Asim(pn  P 13  P22  P24  P33  P44)  (0  +  B  sim^  (216) 

where 

r2fu  2f13  0  0  0  0  > 

f31  f|l+f33  0  0  fl3  0 

=  0  0  2f22  2f24  0  0  B 

sim“  0  .0  f42  f22  +  f44  0  f24  ’ 

0  2f31  0  0  2f33  0 

^  0  0  0  2f42  0  2  f„j 

and  u  =  1  V  time.  The  linear  system  in  eq.  (216)  can  then  be  simulated  which  results  in 
the  time  history  of  the  covariance  matrix.  The  MATLAB  command  LSIM  [15]  can  be 
used,  along  with  other  routines,  to  solve  eq.  (216). 

Having  obtained  the  covariance  matrix  for  each  time  instant,  a  phase  screen  [23] 
must  be  generated.  This  screen  represents  the  random  effects  of  the  atmosphere  on  a 
wavefront  and  provides  random  draws  of  wavefront  phase  for  each  time  instant.  It  is 
desired  that  the  phase  screens  have  the  correct  spatial-temporal  statistics;  that  is,  those 
given  by  the  covariance  matrices.  To  achieve  this,  consider  the  Cholesky  factorization 
[16]  ofP(t), 

P(t)  =  Rp(t)R;(t),  (218) 

which  can  be  performed  since  the  covariance  matrix  is  guaranteed  to  be  positive  definite. 
Consider  now  a  random  draw  of  a  vector  of  zero-mean,  uncorrelated,  unit-variance, 
Gaussian  random  variables,  b(t).  The  covariance  of  b(t)  is 

E{b(t)bT(t)}  =  I  (219) 
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where  I  is  the  identity  matrix. 

Recall  that  a  realization  of  the  wavefront  phase  can  be  written  as 
0(Rx,  t)  =  ^  a.t  (t)Zj  (x) .  Now  consider  a  random  draw  of  the  vector 

a(t)  =  [a2(t)  a3 (t)  •••  an+1(t)] .  Note  that  ai(t)  is  not  included  since  only  piston 

removed  quantities  are  considered.  This  random  draw  is  accomplished  using  the  vector 
b(t)  in  the  following  manner: 

a(t)  =  Rp  (t)b(t) .  (220) 

Then,  the  covariance  of  a(t)  will  match  the  desired  covariance  given  by  P(t).  This  can  be 
seen  by  considering 

E{a(t)aT(t)}  =  E{Rp(t)b(t)(Rp(t)b(t))T}  =  E{Rp(t)b(t)bT(t)R^(t)} 

=  Rp(t) E{b(t)bT(t)}Rp (t)  =  Rp (t) IRp (t)  =  Rp (t)Rp (t)  =  P(t) .  (221) 

By  repeated  application  of  this  algorithm,  a  phase  screen  is  generated  which  has  the 
desired  covariance  time  history.  The  vector  a(t),  for  each  time  instant,  provides  the  truth 
model  for  the  AO  system.  In  the  development  of  the  phase  screen,  the  covariance  matrix, 
not  the  covariance  kernel  matrix,  was  used.  Recall  that  in  deriving  the  atmospheric 
model  in  Chapter  4,  the  covariance  kernel  matrix  was  used.  Therefore,  the  atmospheric 
model  contains  the  kernel  information.  For  the  phase  screen,  the  Cholesky  factorization 
requires  a  positive  definite  matrix.  Since  the  covariance  matrix  is  positive  definite,  it  was 
decided  to  use  this  information  instead  of  the  covariance  kernel  matrix;  however,  the 
covariance  kernel  matrix  information  is  embedded  in  the  phase  screen  since  this  was  used 
to  determine  the  atmospheric  model. 
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Outputs  from  the  simulation  consist  of  time  histories  of  the  state  variables,  the 
filter’s  estimate  of  these  states,  and  filter  covariance  matrices.  Post-simulation  data 
processing  produces  plots  of  the  outbound  wavefront’s  RMS  phase  distortion,  the  RMS 
phase  distortion  after  wavefront  correction,  and  the  Strehl  ratios.  In  addition  to  these 
single  realizations,  a  Monte  Carlo  analysis  consisting  of  10  runs  was  conducted  to 
generate  ensemble  statistics.  In  terms  of  ensemble  statistics,  the  following  equations 
were  used: 


where 


and 


i  NR 

Xek(0’  (222) 

JNK  k=l 

1  NR 

s«(t)  =  —  -  X[ek(t)-mc(t)][ek(t)-me(t)]  ,  (223) 

JNK  -  1  k=l 

NR  =  number  of  Monte  Carlo  runs, 
ek(t)  =  ak(t)-ak(t)  =  error, 
k  =  sample  number  realization, 
me(t)  =  mean  of  random  process  ek(t) , 

Sg  (t)  =  variance  of  random  process  ek  (t) . 


7.4  Simulation  Results 

In  this  section,  the  plots  obtained  from  the  simulation  will  be  displayed.  Since  it 
is  expected  that  shot  noise  and  photon  count  entering  the  HWFS  will  vary  greatly  and 
these  are  the  main  components  of  measurement  noise,  a  parameter  study  was  performed 
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on  the  noise  variance.  The  robustness  of  the  controller,  to  a  parameter  variation,  was  also 
evaluated  by  varying  the  process  noise  in  the  truth  model.  Recall,  when  synthesizing  the 
atmospheric  model,  that  the  process  noise  strength  was  set  to  Qa  =  I.  In  the  simulations, 
Qa  of  the  truth  model  was  set  to  I,  101,  and  201.  For  the  larger  strengths,  the  simulation  is 
being  conducted  with  a  mismatched  model.  If  the  system  performs  well  in  these  cases,  it 
can  be  concluded  that  the  complete  model  is  accurate  and  robust  to  variations  in  the  Qa 
parameter.  Lastly,  the  actuator  time  constant  was  varied.  It  is  expected  that  the  actuators 
for  the  ABL  will  be  very  fast.  In  fact,  it  may  be  possible  to  model  the  actuators  as  a 
constant  gain  multiplier,  i.e.,  able  to  respond  instantaneously  to  changes  in  command 
inputs.  However,  in  order  to  simulate  a  basic  time  delay,  the  actuator  time  constant  was 
made  larger  by  five  orders  of  magnitude.  By  no  means  does  this  represent  a  transport 
delay.  In  fact,  the  larger  time  constant  simply  delays  the  application  of  the  conjugate 
phase  to  the  mirror.  Table  3  lists  the  parameters  for  the  simulations. 


Ri 

Qa 

Ta  (sec.) 

SIMULATION  #1 

54.74 

I 

0.0000075 

SIMULATION  #2 

25.16 

I 

0.0000075 

SIMULATION  #3 

25.16 

101 

0.0000075 

SIMULATION  #4 

201 

0.0000075 

SIMULATION  #5 

25.16 

I 

0.75 

Table  3.  Simulation  parameters. 
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In  Table  3,  the  values  of  Rj  were  determined  by  the  MATLAB  [15]  simulation.  The 
desired  values  were  50  and  25.  However,  with  a  finite  number  of  samples,  precise  values 
are  unobtainable. 

7.4.1  Simulation  Run  #1 

The  first  plots  shown  here  are  a  sample  realization  of  the  outbound  wavefront’s  x- 
and  y-tilts  from  the  phase  screen  with  the  filter’s  estimates  superimposed.  In  this  figure, 
the  phase  screens  are  in  red  while  the  filter  estimates  are  in  blue. 

Filter  error  is  defined  as  the  true  state  minus  the  filter’s  estimate  of  the  state  (state 
refers  to  the  Zernike  coefficients).  The  filter  covariance,  Pf,  is  the  filter’s  indication  of 
the  uncertainty  in  its  estimates,  as  in  Eq.  (174).  The  square  root  of  the  (j,  j)  element  of 
the  Pf  matrix  is  what  the  filter  believes  to  be  the  one  sigma  value  of  its  error. 


Outbound  Wavefrontx-tilt:Estimate  &  Screen 


Time  (sec) 

Figure  26.  Outbound  wavefront  tilts:  estimates  and  screen. 


109 


Although  Figure  26  shows  that  the  filter  appears  to  be  tracking  the  true  state,  it  does  not 
indicate  how  the  filter  error  compares  with  what  it  thinks  are  its  one  sigma  values. 
Figures  27  and  28  show  a  single  sample  time  history  of  the  filter  error  for  the  x-  and  y- 
tilts,  along  with  the  filter  computed  one  sigma  values. 


F  i Ite  r  e  rro  r  a  n d  +  -  filte  r  s  tnd  .  d  e  v.,  o  utb  o  u  nd  x-tilt 
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Figure  27.  Outbound  x-tilt  filter  error  and  filter  standard  deviation. 
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Figure  28.  Outbound  y-tilt  filter  error  and  filter  standard  deviation. 
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Looking  at  these  single  sample  realizations  of  filter  error,  it  appears  that 
±^/Pf  (3,3)  and  ±JP(  (4,4)  are  reasonable  one  sigma  values  for  the  eiTor  and  that  the 

error  is  zero-mean.  Note  that  the  (3,3)  and  (4,4)  elements  of  the  Pt  matrix  are  used  since 
the  outbound  wavefront  tilts  are  the  third  and  fourth  components  of  the  state  vector  a(t). 

The  actual  values  for  the  filter  standard  deviations  are  ^Pf~(3,3)  =  0.074  and  ^/pr+ (3,3) 

=  0.020  where  the  superscripts  and  V  designate  values  before  measurement  update  (- 
)  and  after  measurement  update  (+).  Also,  -^Pf”(4,4)  =  0.063  while  -^Pf+(4,4)  =  0.012. 

Time-averaged  means  and  standard  deviations  of  the  filter  errors  were  also  computed. 

For  the  x-tilt  case,  the  mean  error  is  -0.00023,  while  for  the  y-tilt,  the  mean  error  is 
0.00021 .  The  standard  deviation  of  the  x-tilt  error  is  0.045  and  that  of  the  y-tilt  error  is 
0.026.  These  statistics  were  computed  using  the  single  sample  run  errors.  Thus,  they  do 
not  accurately  portray  the  average  system  performance.  In  order  to  get  a  true  description 
of  the  errors,  a  Monte  Carlo  analysis  must  be  conducted. 

In  order  to  obtain  the  true  error  mean  and  variance,  a  Monte  Carlo  analysis 
consisting  of  10  runs  was  performed.  Statistics  were  computed  using  eqs.  (222)  and 
(223).  Figures  29  and  30  show  the  mean  and  mean  ±  one  standard  deviation  of  the  filter 
error  for  the  outbound  wavefront  x-  and  y-  tilts.  It  appears  that  the  error  process  is 
approximately  zero-mean  with  standard  deviations  less  than  the  ±^P,  (3,3)  and 

±^/Pr  (4,4)  values  from  Figures  27  and  28.  Hence,  it  appears  that  the  filter  is  not 

properly  tuned.  The  time-average  of  the  mean  x-tilt  error  is  -0.001 1,  the  time-average  of 
the  mean  y-tilt  error  is  0.00093,  and  the  time-averaged  standard  deviations  are  0.037  for 


ill 


the  mean  x-tilt  error  and  0.019  for  the  mean  y-tilt  error.  Hence,  the  standard  deviations 


are  much  less  than  the  average  -y/pr+(j,  j)  and  ^P,  (j.j)  values.  Thus,  the  filter  is  not 
properly  tuned. 


Outbound  x-tilt  error;  mean  and  mean  +  -one  stnd.deviation 


Figure  29.  Outbound  x-tilt,  mean  and  mean  ±  one  standard  deviation. 


Outbound  y-tilt  error:  mean  and  mean  +  ■  one  stnd.deviation 


Figure  30.  Outbound  y-tilt,  mean  and  mean  ±  one  standard  deviation. 


Figures  31-33  show  single  sample  realizations  of  the  controlled  variables.  Recall 
that  it  is  desired  to  drive  the  controlled  variables  to  zero,  that  is,  to  have  -a2(l,)(t)  =  a2(m)(t) 
and  -a3(o)(t)  =  a3(m)(t).  Notice  that  there  is  a  little  error  in  the  x-tilt  case,  but  almost  zero 
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error  in  the  y-tilt  case.  This  is  because  of  the  planar  engagement  geometry  discussed  in 
Chapter  3.  In  this  engagement,  dynamic  motion  occurs  only  in  the  x-direction.  Figures 
32  and  33  show  a  sample  realization  of  the  controlled  variables  before  (blue)  and  after 
(red)  correction. 


Negative  M  irror  and  0  utbound  W  avefront  x-tilt 


Figure  31 .  Controlled  variables. 


x-tiltcontrolled  variables  withand  withoutcompensationvs.time 


Figure  32.  X-tilt  controlled  variables. 
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y  -  tilt  controlled  variables  with  and  without  compensation  vs.  time 


Time  (sec) 

Figure  33.  Y-tilt  controlled  variables. 

It  appears  that,  after  compensation,  the  tilt  variables  are  much  nearer  their  desired  set 
point,  namely,  zero;  however,  ensemble  statistics  must  also  be  computed.  Figures  34  and 
35  show  the  mean  and  mean  ±  one  standard  deviation  of  the  x-  and  y-tilt  controlled 
variables. 


X  -  tilt  controlled  variables:  mean  and  mean  +  -  one  stnd.deviation 


Time  (sec) 

Figure  34.  X-tilt  controlled  variables:  mean  and  mean  ±  one  standard  deviation. 
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Y  -  tilt  control  led  variables:  mean  and  mean  +  -  one  stnd.  deviation 


Figure  35.  Y-tilt  controlled  variables:  mean  and  mean  ±  one  standard  deviation. 

From  these  figures,  it  can  be  seen  that  the  controlled  variables  are  nearly  zero-mean.  The 
y-tilt  variance  is  smaller  than  that  of  the  x-tilt.  Again,  this  can  be  attributed  to  the  planar 
engagement  geometry.  The  time-averaged  mean  and  standard  deviation  of  the  x-tilt 
controlled  variables  are  -0.00096  and  0.031,  while  the  time-averaged  mean  and  standard 
deviation  of  the  y-tilt  controlled  variables  are  -0.00017  and  0.023.  Thus,  the  standard 
deviation  on  the  y-tilt  controlled  variables  is  less  than  those  of  the  x-tilt  controlled 
variables,  with  both  processes  being  nearly  zero-mean. 

Figures  36  and  37  show  the  Strehl  ratio  after  tilt  compensation  minus  the 
outbound  wavefront  Strehl  ratio  and  the  simple  conjugation  Strehl  ratio,  respectively, 
along  with  a  zero  line  for  reference.  These  figures  were  generated  using  a  single  sample 
realization. 
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Strehl  after  ti It  comp  minus  strehl  before  correction 


Figure  36.  Strehl  ratio  after  tilt  compensation  minus  Strehl  ratio  before  correction. 


Strehl  after  ti  It  comp  minus  strehl  after  simple  conjugation 


Figure  37.  Strehl  ratio  after  tilt  compensation  minus  Strehl  ratio  after  simple 

conjugation. 

Figures  36  and  37  show  the  improvement  in  Strehl  ratio  by  using  the  tilt  compensation 
scheme.  From  the  last  two  figures,  it  can  be  seen  that,  for  this  particular  single  sample 
realization,  the  Strehl  is  increased  by  an  average  of  0.144  for  the  no  compensation  case 
and  about  0.031  for  the  simple  conjugation  case. 
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Figures  38-40  show  the  mean  and  mean  ±  one  standard  deviation  of  the  three 
Strehl  ratios.  The  mean  Strehl  is  0.484  for  the  uncompensated  case,  0.638  for  tilt 
compensation,  and  0.591  for  simple  conjugation.  Hence,  better  Strehl  ratios  are  obtained 
when  the  tilt  compensation  scheme  is  used.  Also,  the  standard  deviation  decreases 
slightly  from  the  uncompensated  case  (standard  deviation  of  0.293),  to  simple 
conjugation  (standard  deviation  =  0.267),  to  the  tilt  compensation  case  (standard 
deviation  =  0.245). 

Strehl  ratio  before  compensation:  mean  and  mean  +  -  o  n  e  stnd.  deviation 


Time  (sec) 


Figure  38.  Uncompensated  Strehl  ratio. 


Strehl  ratio  aftertilt  compensation:  mean  and  mean  +  •  one  stnd.  deviation 


T  Im  e  (sec) 


Figure  39.  Tilt  compensation  Strehl  ratio. 
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Strehl  ratio  after  simple  conjugation:  mean  and  mean  +  -  one  stnd.  deviation 


0  0  .2  0  .4  0  .6  0  .8  1  1  .2 

T  Im  e  (sec) 


Figure  40.  Simple  phase  conjugation  Strehl  ratio. 


Figure  41  shows  a  single  sample  realization  of  the  RMS  phase  distortion  of  the 
outbound  (red)  and  compensated  (blue)  wavefronts. 


Outbound  and  corrected  RMS  phase  distortion 


Time  (sec) 


Figure  41.  RMS  phase  distortions. 


It  appears  that  the  RMS  phase  distortion  is  decreased  using  tilt  compensation.  Figures 
42  and  43  show  the  RMS  statistics  from  the  Monte  Carlo  analysis.  It  can  be  seen  that 
the  mean  and  variance  are  significantly  reduced  by  employing  tilt  compensation. 
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Outbound  Wave  front  RM  S  :  m  ean  and  mean+-one  stnd.  deviation 


Figure  42.  Outbound  wavefront  RMS  statistics. 


Corrected  Wave  front  RMS:  mean  and  mean  +  -  o  n  e  stnd.  deviation 


Figure  43.  Corrected  wavefront  RMS  statistics. 


For  the  outbound  wavefront,  the  time-averaged  mean  and  standard  deviation  are  0.102 
and  0.055,  respectively.  For  the  corrected  wavefront,  the  time-averaged  mean  and 
standard  deviation  are  0.048  and  0.024.  Hence,  the  RMS  means  and  variances  are 
reduced  using  tilt  compensation. 

In  conclusion,  the  system  seems  to  be  performing  adequately  for  this  particular 
set  of  parameters,  i.e.,  measurement  and  process  noise  strength  and  actuator  time 
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constant.  That  is,  the  Strehl  ratios  and  RMS  phase  distortions  are  notably  better  when 
tilt  compensation  is  applied.  Table  4  displays  the  time-averaged  means  and  standard 
deviations  of  the  plots  displayed  in  this  section. 

However,  since  the  measurement  noise  variance  is  unknown,  more  simulations 
must  be  run.  Simulation  set  #2  examines  the  effects  of  decreasing  the  measurement 
noise  variance.  In  this  case,  the  filter  will  weigh  the  measurements  more  heavily  and 
place  less  emphasis  on  the  internal  model. 


Figure  # 

Variable  Name 

Mean 

Std.  Deviation 

27 

X-tilt  filter  error 

-0.00023 

0.045 

28 

Y-tilt  filter  error 

0.00021 

0.026 

29 

Outbound  x-tilt 

-0.0011 

0.037 

30 

Outbound  y-tilt 

0.00093 

0.019 

34 

X-tilt  controlled  variables 

-0.00096 

0.031 

35 

Y-tilt  controlled  variables 

-0.00017 

0.023 

38 

Uncompensated  Strehl  ratio 

0.484 

39 

Tilt  compensation  Strehl  ratio 

0.638 

0.245 

40 

Simple  phase  conjugation  Strehl  ratio 

0.591 

42 

Outbound  wavefront  RMS 

0.102 

| 

43 

Corrected  wavefront  RMS 

0.048 

0.024 

Table  4.  Statistics  for  run  #1  plots. 


7.4.2  Simulation  Run  #2 

In  this  case,  R;  is  nearly  one-half  that  of  run  #1 .  Figure  44  shows  a  single  sample 
realization  of  tilts  while  Figures  45  and  46  show  single  sample  realizations  of  filter  error 
and  filter  computed  one  sigma  values. 
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Outbound  Wavefront  x-tilt:  Estimate  &  Screen 


Figure  44.  Outbound  wavefront  tilts:  estimates  and  screen. 


The  mean  errors  are  0.00019  and  -0.0002  for  the  x-  and  y-tilts,  respectively.  In  this  case, 
the  filter  error  standard  deviations  are  smaller  than  those  of  the  filter  errors  shown  in 
Figures  27  and  28.  Recall  that  in  simulation  run  #1,  the  standard  deviations  were  0.045 
(x-tilt)  and  0.026  (y-tilt).  Here,  the  standard  deviation  of  the  x-tilt  filter  error  is  0.033 

while  that  of  the  y-tilt  is  0.018.  The  filter  computed  standard  deviations  are:  ^/pf+(3,3)  = 

0.014,  ^Pf~(3,3)  =0.052,  ^/pf+(4.4)  =  0.008,  and  ^Pr'f44)  =0.045.  As  compared  to 

simulation  run  #1,  the  filter  computed  standard  deviations  are  less  which  implies  that  the 
filter  is  more  confident  about  its  estimates. 


121 


Filter  error  and  +  -  filter  stnd  dev.,  outbound  x-tilt 
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Figure  45.  Outbound  x-tilt  filter  error  and  filter  standard  deviation. 
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Filter  error  and  +  -  filte  r  stnd.  dev.,  outbound  y-tilt 


Time  (sec) 


Figure  46.  Outbound  y-tilt  filter  error  and  filter  standard  deviation. 


Figures  47  and  48  show  the  mean  and  mean  ±  one  standard  deviation  of  the 
outbound  tilts.  The  errors  are  approximately  zero-mean  and  once  again,  the  variances  are 
smaller  than  those  displayed  in  simulation  run  #1 .  For  run  #2,  the  standard  deviations  of 
the  error  are  0.03  (x-tilt)  and  0.017  (y-tilt).  Those  for  run  #1  were  0.037  (x-tilt)  and 

0.019  (y-tilt).  Also,  the  standard  deviations  are  nearly  equal  to  the  P,  (3,3)  and 
±^/Pf  (4,4)  values.  This  would  imply  that  the  filter  is  better  tuned. 
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Outbound  x-tilt  e  rro  r:  mean  and  mean  +- one  stnd.  deviation 


Figure  47.  Outbound  x-tilt,  mean  and  mean  ±  one  standard  deviation. 


Outbound  y-tilt  error:  m  ean  and  m  ean  +  -  one  stnd.  deviation 


Figure  48.  Outbound  y-tilt,  mean  and  mean  ±  one  standard  deviation. 


Figure  49  shows  a  single  sample  realization  of  the  controlled  variables,  while 
Figures  50  and  5 1  show  the  mean  and  mean  ±  one  standard  deviation  of  the  controlled 
variables.  Here,  the  means  and  variances  are  smaller  than  those  displayed  in  run  #1 .  The 
time-averaged  means  of  the  controlled  variables  are  as  follows  (those  from  run  #1  are  in 
parenthesis):  -0.00016  (-0.00096)  for  x-tilt  and  0.00015  (-0.00017)  for  y-tilt.  The  time- 
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m-x 


averaged  standard  deviations  are  (run  #1  values  in  parenthesis):  0.027  (0.031)  for  x-tilt 
and  0.017  (0.023)  for  y-tilt. 


Negative  Mirror  and  Outbound  Wavefrontx-tilt 


Figure  49.  Controlled  variables. 


X  -  tilt  controlled  variables  vs.  tim  e:  m  ean  and  mean  +  ■  one  stnd.  deviation 


Figure  50.  X-tilt  controlled  variables. 


124 


Y  -  tilt  controlled  variables  vs.  time:  mean  and  mean  +  •  one  stnd.  deviation 


Figure  51 .  Y-tilt  controlled  variables. 


Figures  52  and  53  show  the  differences  in  Strehl  ratios  from  a  single  sample 
realization.  There  is  an  increase  in  the  gain  in  Strehl  in  this  situation  as  compared  to  run 
#1 .  Hence,  the  controller  seems  to  perform  better  with  this  set  of  simulation  parameters. 
This  performance  should  be  expected  since  the  measurements  are  more  precise.  In  fact, 
the  average  change  in  Strehl  ratio  is  0.164  for  no  compensation  and  0.038  for  simple 
conjugation.  For  run  #1,  the  increases  were  0.144  and  0.031. 


Strehl  after  tilt  comp  minus  strehl  before  correction 


Figure  52.  Strehl  ratio  after  tilt  compensation  minus  Strehl  ratio  before  correction. 
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Strehl  after  tilt  comp,  minus  strehl  after  simple  conjugation 


Figure  53.  Strehl  ratio  after  tilt  compensation  minus  Strehl  ratio  after  simple 

conjugation. 


Figures  54-56  show  the  mean  and  mean  ±  one  standard  deviation  of  the  Strehl  ratios. 


Strehl  ratio  before  compensation:  mean  and  mean  +  -  o  n  e  stnd.  deviation 


Figure  54.  Uncompensated  Strehl  ratio. 


Strehl  ratio  afte  r  tilt  com  p  ensation:  mean  and  mean  +  ■  one  stnd.  deviation 


Figure  55.  Tilt  compensation  Strehl  ratio. 
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Strehlratio  after  simple  conjugation:  mean  and  mean+-one  stnd.  deviation 
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Figure  56.  Simple  phase  conjugation  Strehl  ratio. 
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The  time-averaged  mean  Strehl  ratios  are  (run  #1  values  in  parenthesis):  0.481 
uncompensated,  0.621  (0.591)  simple  phase  conjugation,  and  0.643  (0.638)  for  tilt 
compensation.  Hence,  the  Strehl  ratios  have  increased.  The  time-averaged  standard 
deviations  are  (run  #1  values  in  parenthesis)  are:  0.285  uncompensated,  0.265  (0.267) 
simple  conjugation,  and  0.242  (0.245)  for  tilt  compensation.  Thus,  the  Strehl  ratios  have 
increased  and  the  variances  decreased  as  compared  to  run  #1 . 

Figures  57-59  show  the  RMS  phase  distortion  and  the  associated  Monte  Carlo 
statistics.  Note  that  the  corrected  RMS  values  are  smaller  than  those  in  run  #1 .  Recall 
that  with  Rj  =  25.16  instead  of  54.74,  more  emphasis  is  placed  on  the  measurements  and 
less  is  placed  on  the  internal  model. 
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Corrected  RMS  Outbound  RMS 


0.3  5 


Outbound  and  corrected  RMS  phase  distortion 


Figure  57.  RMS  phase  distortions. 


Outbound  Wave  front  R  M  S  :  m  ean  and  mean  +  -  one  stnd.  deviation 


Figure  58.  Outbound  wavefront  RMS  statistics. 


Corrected  Wave  front  RMS:  mean  and  mean  +■  one  stnd.  deviation 


Figure  59.  Corrected  wavefront  RMS  statistics 


The  time-averaged  mean  RMS  values  are  (run  #1  values  in  parenthesis):  0.104  for  the 
outbound  wavefront  and  0.045  (0.048)  for  the  corrected  wavefront.  The  time-averaged 
RMS  standard  deviations  are:  0.051  (0.055)  for  the  outbound  wavefront  and  0.021 
(0.024)  for  the  corrected  wavefront. 

In  this  case,  system  performance  is  better  than  in  run  #1 .  That  is,  more  accurate 
estimates  of  the  outbound  wavefront  Zemike  coefficients  are  made.  This,  in  turn,  results 
in  better  regulation  of  the  controlled  variables  and  hence,  better  Strehl  ratios  and  RMS 
phase  distortions.  Table  5  summarizes  the  results  of  this  simulation  run. 


IBISWRI-l 

Variable  Name 

Mean 

Std.  Deviation 

45 

X-tilt  filter  error 

0.00019 

0.033 

46 

Y-tilt  filter  error 

-0.0002 

0.018 

47 

Outbound  x-tilt 

-0.00022 

0.03 

48 

Outbound  y-tilt 

0.00015 

0.017 

X-tilt  controlled  variables 

-0.00016 

0.027 

51 

Y-tilt  controlled  variables 

0.00015 

0.017 

54 

Uncompensated  Strehl  ratio 

0.481 

0.285 

55 

Tilt  compensation  Strehl  ratio 

0.643 

0.242 

56 

Simple  phase  conjugation  Strehl  ratio 

0.621 

0.265 

58 

Outbound  wavefront  RMS 

0.104 

0.051 

59 

Corrected  wavefront  RMS 

0.045 

0.021 

Table  5.  Statistics  for  run  #2  plots. 


Now,  the  robustness  of  the  Kalman  filter  to  variations  in  the  process  noise  strength  will 
be  tested.  This  simulation  is  used  to  determine  the  effects  of  model  mismatches  upon  the 
operation  of  the  system.  Simulation  run  #3  displays  these  results. 
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7.4.3  Simulation  Run  #3 


Recall  that  in  this  case,  the  process  noise  of  the  truth  model,  Qa,  was  set  to  101, 
while  the  filter  was  run  with  Qa  =  I.  Therefore,  the  simulation  is  being  run  to  evaluate 
the  robustness  of  the  filter  to  variations  in  the  process  noise  strength.  Figure  60  shows  a 
sample  realization  of  the  outbound  tilts. 
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Outbound  Wavefrontx-tilt:  Estimate  &  Screen 
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Figure  60.  Outbound  wave  front  tilts:  estimates  and  screen. 


It  still  appears  that  the  filter  is  accurately  estimating  the  outbound  wavefront  coefficients. 
Figures  61  and  62  show  a  single  realization  of  filter  error  along  with  the  filter  computed 
one  sigma  values.  The  mean  errors  are  -0.00021  for  x-tilt  and  0.00041  for  y-tilt.  Thus, 
the  errors  are  still  nearly  zero-mean;  however,  the  standard  deviations  of  the  errors  are 
much  greater  than  in  the  previous  simulation  run.  In  fact,  the  x-tilt  standard  deviation  is 
0.057  while  the  y-tilt  standard  deviation  is  0.036.  For  run  #2,  the  standard  deviations 
were  0.033  (x-tilt)  and  0.01 8  (y-tilt).  Hence,  the  model  inadequacies  are  causing 
problems  with  the  estimation  of  the  outbound  wavefront  Zemike  coefficients.  In  other 
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words,  the  filter  underestimates  the  size  of  its  own  errors  since  Qa  =  10 1  for  the  truth 


model  and  Qa  =  I  for  the  filter. 


Filter  error  and  +  -  filter  stnd  .  dev.,  outbound  x-tilt 


Time  (sec) 


Figure  61 .  Outbound  x-tilt  filter  error  and  filter  standard  deviation. 


Filter  error  and  +  -  filter  stnd.  dev.,  outbound  y-tilt 
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Figure  62.  Outbound  y-tilt  filter  error  and  filter  standard  deviation. 


Figures  63  and  64  show  the  mean  and  mean  ±  one  standard  deviation  of  the  error  in  the 
outbound  tilts.  The  errors  are  zero-mean,  but  the  error  variances  are  larger  than  those  in 
run  #2.  The  error  standard  deviations  for  this  run  are  (run  #2  values  in  parenthesis): 
0.054  (0.03)  for  the  x-tilt  and  0.032  (0.017)  for  the  y-tilt. 
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Outbound  x-tilterror:mean  and  mean  +  ■  one  stnd.  deviation 


Time  (sec) 


Figure  63.  Outbound  x-tilt  error. 


Outbound  y-tilterror:mean  and  mean  +  -  one  stnd.deviation 


Time  (sec) 


Figure  64.  Outbound  y-tilt  error. 


Figures  65-67  display  the  controlled  variables.  The  controlled  variables  are 
approximately  zero-mean,  however,  the  variances  are  significantly  greater  as  compared 
to  run  #2. 
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Negative  Mirror  and  OutboundWavefrontx-tilt 


Figure  65.  Controlled  variables. 


X  -  tilt  controlled  variables  vs.  tim  e:  m  ean  and  mean  +  •  one  stnd.deviation 


Figure  66.  X-tilt  controlled  variables. 


Y  -  ti  It  control  led  variables  vs.  time:  mean  and  mean  +  -  one  stnd.deviation 


Figure  67.  Y-tilt  controlled  variables. 
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Although  the  standard  deviations  of  the  controlled  variables  are  larger  than  those 
displayed  in  simulation  run  #2,  the  controller  is  still  regulating  the  controlled  variables  to 
zero.  The  standard  deviations  of  the  controlled  variables  are  (run  #2  values  in 
parenthesis):  0.039  (0.027)  for  the  x-tilt  and  0.031  (0.017)  for  the  y-tilt. 

Figures  68-70  show  ensemble  statistics  of  the  Strehl  ratios.  Although  the 
compensated  Strehl  ratios  are  still  high,  they  are  worse  than  in  the  previous  two  cases. 
The  tilt  compensated  time-averaged  mean  Strehl  ratio  is  0.495,  while  for  run  #2  the  value 
was  0.643.  The  simple  phase  conjugation  Strehl  ratio  is  0.47,  while  for  run  #2  it  was 
0.621 .  However,  the  compensated  Strehls  are  better  than  the  uncompensated  case  (Strehl 
of  0.44).  This  would  imply  that  even  with  model  inadequacies  of  the  magnitude 
described  in  this  simulation  run,  tilt  compensation  and  simple  conjugation  are  better  than 
no  compensation. 
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Strehl  ratio  before  compensation:  mean  and  mean+-one  stnd.  deviation 


Figure  68.  Uncompensated  Strehl  ratio. 
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Strehl  ratio  after  tilt  compensation:  mean  and  mean  +  -  one  stnd.  deviation 


Figure  69.  Tilt  compensation  Strehl  ratio. 


Figure  70.  Simple  conjugation  Strehl  ratio. 


System  performance  is  definitely  worse  in  this  case  as  compared  to  run  #2, 
however,  it  appears  that  the  controller  is  rather  robust  to  variations  in  the  Qa  parameter. 
That  is,  the  Strehl  ratios  are  higher  for  the  tilt  compensation  and  simple  conjugation 
cases  as  compared  to  no  compensation.  Hence,  even  if  the  atmospheric  model  does  not 
provide  a  precise  representation  of  the  physical  plant,  system  performance  can  be 
enhanced  by  applying  tilt  compensation.  Table  6  summarizes  this  simulation  run. 


135 


Figure  # 

Variable  Name 

Mean 

Std.  Deviation 

61 

X-tilt  filter  error 

-0.00021 

62 

Y-tilt  fdter  error 

0.00041 

0.036 

63 

Outbound  x-tilt 

-0.00029 

64 

Outbound  y-tilt 

0.00026 

0.032 

66 

X-tilt  controlled  variables 

-0.00082 

67 

Y-tilt  controlled  variables 

-0.00021 

0.031 

68 

Uncompensated  Strehl  ratio 

0.44 

— 

69 

Tilt  compensation  Strehl  ratio 

0.495 

— 

70 

Simple  phase  conjugation  Strehl  ratio 

0.47 

— 

Table  6.  Statistics  for  run  #3  plots. 


7.4.4  Simulation  Run  #4 

In  this  case,  Qa  of  the  truth  model  was  set  to  201.  Hence,  the  model  uncertainties 
are  more  severe  than  those  in  simulation  run  #3.  Figure  71  shows  a  single  sample 
realization  of  the  x-  and  y-  tilts. 
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Figure  7 1 .  Outbound  wavefront  tilts:  estimate  and  screen. 
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Figures  72  and  73  show  a  single  realization  of  filter  error  along  with  the  filter  computed 
one  sigma  values  while  Figures  74  and  75  show  the  mean  and  mean  +  one  standard 
deviation  of  the  error  in  the  outbound  tilts. 


Filter  error  and  +  -  filter  stnd.  dev.,  outbound  x - 1 i ft 


Figure  72.  Outbound  x-tilt  filter  error  and  filter  standard  deviation. 


Filter  error  and  +  -  filter  stnd.  dev.,  outbound  y-tilt 


Figure  73.  Outbound  y-tilt  filter  error  and  filter  standard  deviation. 


The  mean  x-tilt  filter  error  is  -0.0003,  while  the  mean  y-tilt  filter  error  is  0.00013.  The 
standard  deviation  of  the  x-tilt  filter  error  is  0.07  and  the  standard  deviation  of  the  y-tilt 
filter  error  is  0.045.  Recall  that  these  standard  deviations  were  0.033  (run  #2  x-tilt). 
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0.057  (run  #3  x-tilt),  0.018  (run  #2  y-tilt),  and  0.036  (run  #3  y-tilt).  Hence,  the  standard 
deviations  are  larger  which  implies  more  uncertainty  in  the  errors. 


Outbound  x-tilt  error:  mean  and  mean  +-  one  stnd.deviation 
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Figure  74.  Outbound  x-tilt  error. 


Outbound  y-tilt  error:  mean  and  mean  +  -  one  stnd.deviation 


Figure  75.  Outbound  y-tilt  error. 


Also,  the  outbound  wavefront  tilt  errors  are  much  larger  than  in  the  previous  cases.  The 
time-averaged  standard  deviation  of  the  mean  x-tilt  error  is  0.067  and  the  time-averaged 
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standard  deviation  of  the  mean  y-tilt  error  is  0.038.  Previously,  these  values  were  0.03 
(run  #2  x-tilt),  0.054  (run  #3  x-tilt),  0.017  (run  #2  y-tilt),  and  0.032  (run  #3  y-tilt). 
Hence,  the  model  inaccuracies  are  degrading  the  performance  of  the  controller;  however, 
the  errors  are  still  zero-mean. 

Figures  76  and  77  show  the  mean  and  mean  ±  one  standard  deviation  of  the 
controlled  variables.  Although  the  mean  and  variances  have  increased  from  the  case 
when  Qa  =  101,  the  controlled  variables  are  zero-mean. 


X  -ti  It  controlled  variables  vs.  time:  mean  and  mean  +  ■  one  stnd.  deviation 
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Figure  76.  X-tilt  controlled  variables. 


Y  -  tilt  controlled  variables  vs.  time:  mean  and  mean  +  -  one  stnd.deviation 


0  .6 

Time  (sec) 


Figure  77.  Y-tilt  controlled  variables. 
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The  standard  deviation  of  the  x-tilt  controlled  variables  is  0.064  and  for  the  y-tilt  the 
standard  deviation  is  0.04.  For  the  last  two  runs,  the  standard  deviations  were  0.027  (run 
#2  x-tilt),  0.017  (run  #2  y-tilt),  0.039  (run  #3  x-tilt),  and  0.031  (run  #3  y-tilt). 

Figures  78-80  show  the  three  Strehl  ratios.  In  this  case,  the  compensated  Strehl 
values  are  about  the  same  as  the  uncompensated  Strehl  values.  The  time-averaged  mean 
uncompensated  Strehl  ratio  is  0.43 1,  while  the  time-averaged  means  of  the  tilt 
compensation  and  simple  conjugation  Strehl  ratios  are  0.444  and  0.421,  respectively. 
From  run  #3,  the  Strehl  ratios  were  0.495  (tilt  compensated)  and  0.47  (simple 
conjugation). 

Strehl  ratio  before  compensation:  mean  and  mean+-one  stnd.  deviation 


Figure  78.  Uncompensated  Strehl  ratio. 
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S  trehl  ratio  after  tilt  com  pensation:  mean  and  mean  +  ■  one  stnd.deviation 


Time  (sec) 

Figure  79.  Tilt  compensation  Strehl  ratio. 


Strehlratio  aftersimple  conjugation:  mean  and  mean +- one  stnd.deviation 


0  0.2  0.4  0.6  0.8  1  1.2 

Tim  e  (sec) 


Figure  80.  Simple  conjugation  Strehl  ratio. 

In  conclusion,  the  severe  model  inaccuracies  investigated  in  this  simulation  run 
do  hamper  the  system’s  performance.  In  a  number  of  cases,  the  filter’s  errors  are  large. 
These  have  the  effect  of  applying  the  incorrect  phase  compensation  to  the  deformable 
mirror.  Hence,  at  some  time  instants,  the  Strehl  ratios  can  actually  be  less  when  applying 
compensation  with  these  model  inadequacies.  Table  7  displays  the  results  of  this 
simulation  run. 
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Figure  # 

Variable  Name 

Mean 

Std.  Deviation 

72 

X-tilt  filter  error 

-0.0003 

0.07 

73 

Y-tilt  filter  error 

0.00013 

0.045 

74 

Outbound  x-tilt 

0.00026 

0.067 

75 

Outbound  y-tilt 

-0.00017 

0.038 

76 

X-tilt  controlled  variables 

-0.00071 

0.064 

77 

Y-tilt  controlled  variables 

-0.00029 

0.04 

78 

Uncompensated  Strehl  ratio 

0.431 

— 

79 

Tilt  compensation  Strehl  ratio 

— 

80 

Simple  phase  conjugation  Strehl  ratio 

0.421 

— 

Table  7.  Statistics  for  run  #4  plots. 


7.4.5  Simulation  Run  #5 

In  this  scenario,  the  actuator  time  constant  was  set  to  0.75  sec.  This  simulates  a 
coarse  approximation  to  a  time  delay.  In  fact,  the  delay  occurs  in  the  actuation  of  the 
DM.  Because  the  time  delay  only  affects  actuation  of  the  DM,  the  only  variables  that 
change  in  this  simulation  run  are  the  controlled  variables,  and  the  ones  they  affect.  That 
is,  the  outbound  wavefront  tilts  do  not  change  since  they  depend  only  on  the  performance 
of  the  filter.  However,  Strehl  ratios  and  controlled  variables  are  affected  since  they 
require  actuation  of  the  DM.  Thus,  only  those  plots  which  change  will  be  shown.  For 
the  others,  see  simulation  run  #2. 

The  first  plot  shown  is  a  single  sample  realization  of  the  controlled  variables. 
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Negative  Mirrorand  Outbound  W  a  ve  fro  nt  x-tilt 


Figure  81.  Controlled  variables. 


It  appears  that  there  is  much  more  error  in  the  controlled  variables  as  compared  to  those 
displayed  in  run  #2.  Figures  82  and  83  display  ensemble  statistics  of  the  controlled 
variables. 


X  -  tilt  controlled  variables  vs.  time:  mean  and  mean  +  -  o  n  e  stnd.  deviation 


Figure  82.  X-tilt  controlled  variables. 


Y  -  tilt  controlled  variables  vs.  time:  mean  and  mean  +  ■  one  stnd.  deviation 


Figure  83.  Y-tilt  controlled  variables. 


Here,  it  can  be  seen  that  the  variances  of  the  controlled  variables  are  much  larger  than  in 
the  previous  cases.  The  time-averaged  standard  deviation  of  the  x-tilt  controlled 
variables  is  0.075  while  for  the  y-tilt  the  time-averaged  standard  deviation  is  0.046. 
Recall  that  in  run  #2,  the  time-averaged  standard  deviations  were  0.027  (x-tilt)  and  0.017 
(y-tilt).  In  fact,  as  will  be  discussed  shortly,  applying  AO  with  slow  actuators  may 
actually  degrade  performance  below  that  which  is  obtained  in  the  uncompensated 
system.  To  evaluate  this  statement,  ensemble  statistics  of  the  Strehl  ratios  will  be 
examined.  Figures  84-86  show  the  mean  and  mean  ±  one  standard  deviation  of  the 
Strehl  ratios. 
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Strehl  ratio  Strehl  ratio  s,rehl  ratio 


Strehlratio  after  tilt  compensation:  mean  and  mean  +  -one  slnd  .  deviation 


Figure  84.  Strehl  ratio  before  compensation. 


Strehlratio  aftertilt  compensation:  mean  and  mean  +  -  one  stnd.deviation 


Figure  85.  Strehl  ratio  after  tilt  compensation. 


Strehlratio  aftersimple  conjugation:  mean  and  mean  +  -  one  stnd.  deviation 


Figure  86.  Simple  conjugation  Strehl  ratio. 
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As  can  be  seen  in  the  last  three  figures,  the  Strehl  ratio,  at  many  time  instants,  is  better  in 
the  uncompensated  case.  The  time-averaged  mean  Strehl  ratios  are  0.44 
(uncompensated),  0.37  (tilt  compensation),  and  0.39  (simple  conjugation).  The  same 
information  is  also  displayed  in  the  RMS  phase  distortion. 

In  conclusion,  with  an  actuator  time  constant  which  is  too  large,  it  appears  that 
using  AO  can  degrade  system  performance.  What  happens  is  as  follows:  measurements 
are  taken  and  the  outbound  wave  front  Zernike  coefficients  are  estimated.  However,  the 
estimated  conjugate  phase  is  not  applied  to  the  DM  until  a  later  time  because  of  the 
actuator  delay.  Hence,  compensation  is  being  calculated  for  one  propagation  path  while 
the  actual  compensation  is  applied  to  a  different  propagation  path.  In  other  words, 
anisoplanatism  is  not  properly  taken  into  account.  Thus,  with  a  large  actuation  delay, 
system  performance  can  actually  be  degraded  by  applying  AO.  Table  8  shows  the 
statistics  of  the  variables  in  this  simulation  run. 


Figure  # 

Variable  Name 

Mean 

Std.  Deviation 

82 

X-tilt  controlled  variables 

-0.00083 

0.075 

83 

Y-tilt  controlled  variables 

-0.00032 

0.046 

84 

Uncompensated  Strehl  ratio 

0.44 

—  - 

85 

Tilt  compensation  Strehl  ratio 

0.37 

— 

86 

Simple  phase  conjugation  Strehl  ratio 

0.39 

— 

Table  8.  Statistics  for  run  #5  plots. 
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7.5  Summary 

In  this  chapter,  the  complete  AO  control  system  was  simulated.  A  parameter 
study  was  performed  on  the  measurement  noise  variance  since  it  is  expected  that  this 
value  will  vary  greatly  in  a  practical  situation.  The  process  noise  strength  of  the  truth 
model  and  the  actuator  time  constant  were  also  varied.  In  all  cases,  excluding  the  larger 
time  constant,  the  controller  performs  adequately.  Hence,  accurate  estimates  of  the 
outbound  Zernike  coefficients  are  made  and  the  controlled  variables  are  regulated  near 
zero.  With  a  time  delay  that  is  too  large,  performance  can  be  substantially  degraded. 
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VIII.  Conclusions  and  Recommendations 


8.1  Introduction 

The  principle  objective  in  this  research  was  to  establish  a  unified  approach  to 
beam  control  for  the  application  of  AO  to  the  ABL.  The  objective  was  motivated  by  the 
inadequacy  of  current  analysis  methods  to  address  anisoplanatism,  the  specific  ABL 
engagement  geometry,  and  the  correct  development  of  the  correlation  between  inbound 
and  outbound  wavefronts.  Specifically,  this  research  was  intended  to  assess  the  effects  of 
anisoplanatism,  and  to  quantify  the  performance  improvements  that  can  be  attained  by 
applying  tilt  compensation  to  the  AO  system. 

In  Chapter  1  of  this  dissertation,  an  overview  of  AO  systems  and  the  ABL  was 
presented,  along  with  a  historical  background,  problem  statement,  and  key  results. 
Chapter  2  provided  an  overview  of  atmospheric  models  and  Zernike  polynomials.  The 
ABL  specific  engagement  geometry  was  delineated  in  Chapter  3,  including  all  of  the 
necessary  kinematic  variable  and  geometrical  relationships.  The  correlation  between 
inbound  and  outbound  wavefronts  was  determined  in  Chapter  4,  along  with  calculations 
and  stochastic  modeling.  Chapter  5  presented  the  deformable  and  complete  system 
models,  output  equations,  and  the  equivalent  discrete-time  system.  A  Kalman  filter  and 
LQG  controller  were  constructed  in  Chapter  6,  while  simulations  were  performed  in 
Chapter  7,  including  a  discussion  of  the  simulation  results.  In  the  present  chapter,  the 
significant  advances  of  this  research  are  summarized  and  suggestions  for  future  research 
are  offered. 
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8.2  Significant  Advances 

This  research  has  led  to  the  following  significant  advances  regarding  tilt 

compensation  for  the  ABL. 

•  A  unified  approach  to  applying  tilt  compensation  to  the  ABL,  which  accounts  for 
anisoplanatism  and  system  time  delays. 

•  A  fully  developed  engagement  geometry  for  the  ABL  encompassing  target  and 
aperture  motion  as  well  as  wind  effects. 

•  A  modified  frozen  flow  hypothesis  (FFH)  which  takes  into  account  translation  and 
rotation  of  the  aperture  and  target. 

•  A  new  approach  to  developing  an  atmospheric  model  which  directly  uses  correlation 
kernel  data  instead  of  designing  shaping  filters  for  the  atmospheric  states.  This  results 
in  a  lower-order,  and  therefore,  more  practical  atmospheric  model. 

•  Atmospheric  model  developed  provides  a  good  representation  of  the  atmosphere. 

This  is  validated  by  examining  the  performance  of  the  Kalman  filter  based  upon  this 
model  and  also  upon  less  accurate  representations  of  this  model  (varying  the  process 
noise  strength). 

•  Significant  improvement  in  the  reduction  of  wavefront  phase  deformations  can  be 
obtained  by  simply  performing  tilt  correction.  This  result  is  displayed  using 
simulations  of  the  entire  ABL  AO  control  system. 
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8.3  Summary  of  Results 

The  analysis  conducted  in  the  dissertation  has  led  to  a  number  of  key  results.  A 

summary  of  these  results  is  presented  in  this  section. 

•  Tilt  compensation  can  significantly  improve  the  ABL’s  AO  system.  In  terms  of 
Strehl  ratios,  an  improvement  of  0.15  was  achieved  by  using  tilt  compensation  as 
compared  to  the  case  in  which  no  compensation  was  used.  In  comparison  to  the 
simple  phase  conjugation  case,  the  improvement  was  smaller,  approximately  0.04. 
However,  for  a  Strehl  of  0.5,  a  0.04  increase  results  in  a  8%  improvement. 

•  Tilt  compensation  reduces  the  RMS  phase  distortions  by  about  50%.  Therefore,  it  is 
clear  that  tilts  are  a  major  factor  in  wavefront  phase  deformations. 

•  The  required  control  system  bandwidth  for  this  application  is  greater  than  1  kHz. 
This  is  much  greater  than  the  bandwidth  of  current  astronomical  imaging  systems 
employing  AO.  The  bandwidth  requirement  increases  because  the  control  system 
must  be  updated  at  rates  commensurate  with  the  wavefront  phase  variations. 

•  A  Kalman  filter’s  estimates,  with  an  accurate  internal  model  (plant  and  deformable 
mirror),  can  track  the  wavefront  phase  deformations.  Therefore,  the  effects  of 
anisoplanatism  can  be  taken  into  account  and  accurate  estimates  of  the  outbound 
wavefront  Zemike  coefficients  can  be  made. 

•  The  atmospheric  model  designed  in  this  work  is  robust  to  variations  of  the  process 
noise  strength;  that  is,  the  Kalman  filter  can  provide  reasonable  estimates  of  the 
outbound  wavefront  Zernike  coefficients. 
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8.4  Recommendations  for  Future  Research 

The  research  activity  in  this  dissertation  provided  a  description  of  the  design  of  a 
control  system  for  beam  control  in  the  ABL.  While  this  research  answered  some 
important  questions,  there  are  many  more  that  go  as  yet  unanswered.  There  were  several 
considerations  that  were  outside  the  scope  of  this  work.  These  issues  are  outlined  in  the 
following  paragraphs. 

•  In  this  work,  a  planar  engagement  geometry  was  selected.  In  using  this,  dynamic 
motion  occurred  only  in  the  x-direction.  In  the  future,  a  full  three-dimensional 
geometry  should  be  investigated,  including  dynamic  motion  in  both  the  x-  and  y- 
directions. 

•  Although  tilts  represent  the  majority  of  power  in  wavefront  phase  deformations, 
higher-order  modes  can  be  used  to  capture  some  of  the  remaining  power.  Hence,  a 
model,  using  more  than  two  Zernike  modes,  may  be  useful.  However,  it  should  be 
mentioned  once  again  that  more  modes  implies  a  higher-order  Kalman  filter  which 
requires  more  computation  time.  Some  analysis  was  conducted  using  more  Zernike 
modes  (5  modes  to  be  precise),  but  the  improvements  in  performance  did  not 
outweigh  the  additional  computational  loading. 

•  The  atmospheric  model  generated  here  was  for  one  particular  engagement  geometry 
and  one  set  of  atmospheric  parameters.  A  unified  modeling  approach  would  be 
useful.  Perhaps  a  multiple  model  type  of  adaptive  estimator  [17]  could  be  used.  In 
this  case,  multiple  models  would  be  constructed,  one  for  each  entry  in  a  set  of 
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geometrical  and  atmospheric  parameters.  Then,  using  an  appropriate  probability 
function,  the  filter  estimates  would  be  weighted  to  provide  the  final  result. 

•  A  better  understanding  of  the  WFS  noise  model  would  also  be  advantageous. 

Unfortunately,  shot  noise  and  photon  count,  each  of  which  are  engagement-specific, 
represent  the  majority  of  sensor  noise. 
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Appendix  A:  Evaluation  of  the  Integral  Over  0 


In  this  section,  it  is  desired  to  derive  eq.  (100).  Recall,  from  Chapter  4,  eq.  (97), 


that  the  integral  over  0  was  given  as 

|  e  j27t|Ki|l“lcos(0  e“)[cos|(mf  +mj)0  +  gj  +  coSj(mf -mj)0  +  hj 


i  0a+2ft  i_  |  , -  , 

Jj^J  _J_  f  e~  J2it|K,  ||oc|  cos(  e-©n  ) 
5  2 

e„ 


d0  (224) 


where 


and 


h  =  {x{(‘  -  8™.o)[c-n'  -  *]  -  0 -s^Xc-o'  - 1]} 


(225) 


(226) 


By  a  change  of  variables,  a  trigonometric  identity,  and  the  definition  for  a  Bessel  function 
of  the  first  kind  [28],  denoted  by  J,(*),  it  is  stated,  in  eq.  (100),  that 

3(m,  +rrij) 

INTg  =cos{(mf  +mJ)0a+g}7T(-l)  2  J(mf+raj)(2^|l®l) 


+  cos{(mf  -mj)0a  +h}ji(-l)  2  (27t| 


K, 


a 


(227) 


In  order  to  derive  eq.  (227),  first  consider  a  change  of  variables:  let 

(3  =  0  — 0a=>0  =  p  +  0a,dp  =  d0. 


(228) 


Applying  this  change  of  variables  to  eq.  (224)  gives 


2  K 


INTWl6 


1  I  “*  1  I  -  I 

f  - }2rc  iC|  |d|  cosp 


[•]  d(3 


(229) 
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where 


[•]  =  cosj(mf  +  m j  +  (mf  +  m j  )0a  +  g|  +  cos|(mf  -  rrij  )(3  +  (mf  -  nij  )9a  +  h}  .  (230) 

Considering  the  first  cosine  term  only,  using  the  identity 

cos(a+b)  =  cos(a)cos(b)-sin(a)sin(b),  (231) 

and  identifying  a  =  (mf  +  mj)P,  b  =  (mf  +  mj)0a  +  g,  the  term  can  be  rewritten  as 

cos{(mf  +mj)p  +  (mf  +  m^  +  g}  =  cos{(mf  +mj)p}cos{(mf  +  mj)0a  +g} 

-sin{(mf  +mj)pjsin{(mf  +mj)0a+gj.  (232) 

Likewise,  the  second  cosine  term  in  eq.  (229)  can  be  written  as 

cos|(mf  -  nij  +  (mf  -  )0a  +  h J  =  cosj(mf  -  nij  )p|  cos|(mf  -  nij  )0a  +  hj 

—  sinj(mf  —  mj)pjsinj(mf  -m.  )0a+h}.  (233) 

Therefore,  eq.  (229)  becomes 

INTS  =  f jVJ2^'||S|'">]d|3  (234) 

where 

[•]  =  cos{(mf  +mj)p}cos{(mf  +mj)0a  +g}-sin{(mf  +mj)p}sin{(mf  +mj)ea  +g} 

+cos|(mf  - mjpjcosj^nif  -mj)0a  +h|-sin|(mf  - mj)p|sin|(mf  -mj)0a  +hj.  (235) 
Hence,  there  are  four  integrals  to  evaluate.  Consider  the  first  integral  which  is 
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INTePan.  =~  je^|K,|Hc0sP[cos{(mf  +mJ(3}cos{(mf  +mj)ea  +g}]d(3 


=  — cos|(mf  +mj)0a  +  g}je  j2T‘HC0SP  Cos|(mf  +  mj)p}  dp.  (236) 

2  A 


From  Gradshteyn  and  Ryzhik  [9],  a  Bessel  function  can  be  expressed  as 


1 


3X/ 


27t(-l)  /2  0 


2n 

- 1  e"jzcosp  cos(Xp)dp  . 


(237) 


Identifying  z  =  2rt  k,  a  and  X  =  (mf  +  mj)  and  using  eq.  (237)  in  eq.  (236)  produces 


INT; 


Bpart  1 


=  cos[(mf  +mj)ea+g|7r(-l)  /z  J(m|+mi)(27tK,  |a|j.  (238) 


Using  another  formula  from  Gradshteyn  and  Ryzhik  [9],  namely 

2  K 

J e-jzcosp  sin(^P)dp  =  0,  (239) 

0 

it  can  be  seen  that  the  integrals  with  sine  terms  are  equal  to  zero.  Performing  the  same 
procedure  for  the  remaining  cosine  term  in  eq.  (234)  and  putting  the  results  together  gives 
the  desired  final  form: 


3(ni|  +nij) 

INT§  =  cos{(mf  +mj)ea  +g}7t(-l)  2 

3(m,-mj) 

+cosj(mf  -mj)0a  +h|7t(-l)  2 

which  is  the  same  as  eq.  (100)  in  Chapter  4. 


J 


J 


[mr-mj| 


(240) 
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Appendix  B:  Evaluation  of  Discrete-Time  System,  M  and  N  Matrices,  and  Derivation  of 

RMS  Phase  Distortions 


In  this  section,  it  is  desired  to  evaluate  some  of  the  quantities  discussed  earlier, 
such  as,  the  equivalent  discrete-time  system  matrices  (Section  5.6),  N  matrix  (Section 
5.5),  and  the  M  matrix  (Section  5.2).  Also,  the  details  of  deriving  eq.  (212)  will  be 
shown. 


B.l  Discrete-Time  System  Matrices 

To  begin,  consider  the  equivalent  discrete-time  system  given  by  eq.  (17 1).  In  this 
case,  0(tj+1,  tj),  Bd(tj),  and  Qd(tj)  must  be  delineated.  For  the  simulations,  it  has  been 
determined  that  a  sample  rate  of  tj+i  -  tj  =  0.0008  sec.  is  sufficient.  With  a  constant 
sample  period,  evaluation  of  the  discrete-time  matrices  is  much  simpler.  In  the  case  of 
0(tJ+i,  tj)  and  Qd(tj),  an  algorithm  specified  in  Brown  and  Hwang  [2]  but  developed  by 
van  Loan  [31]  will  be  used.  The  algorithm  is  as  follows: 

1 .  Form  the  matrix 

(v.-‘T  <24» 


-F  GGt 


2. 


Form  eAALG  using  the  MATLAB  command  expm(AALG)  [15],  or  some 
other  equivalent  software  tool: 


Balg  =expm(AALG)  = 


fc-'Qd 

oT 


(242) 
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where  the  time  instants  tj+i  and  tj  have  been  eliminated  since  the  sample 
period  is  a  constant.  The  upper  left  partition  of  BALg  is  of  no  concern. 

3.  Transpose  the  lower  right  partition  of  BAlg  to  obtain  the  state-transition 
matrix. 

4.  Lastly,  Qd  is  obtained  as  follows: 

Qd  =  <I>  *  (upper  right  partition  of  BALg)- 

Bd(tj)  can  be  evaluated  by  solving  the  following  differential  equation  over  each  interval 

[tj,  tj+i)  [16]: 

5(t,tj)  =  FB(t,tj)+B.  (243) 

Equation  (243)  is  integrated  forward  from  B(tj,tj)  =  0  to  time  tj  to  yield 
Bd(tj)  =  B(tj+I ,tj)  [16], 


B.2  Mirror  Matrix  M 

In  this  section,  the  matrix  M  in  eq.  (125)  will  be  developed.  Recall  that  M  is 

M  =  [m2  m3  mn+1]T. 

For  the  case  of  x-  and  y-  tilts,  M  becomes 

M  =  [m2  m3]T  (244) 

where 

m[  =[J\V(x)Zi  (x)I|(x)dx  |w(x)Zj  (x)I2(x)dx  •••  Jw(x)Zj  (x)I„a  (x)dx] .  (245) 
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In  order  to  develop  this  matrix,  the  influence  functions  must  be  specified.  There  are 
many  different  suitable  influence  functions,  i.e.,  pyramids,  cones,  etc.  A  right  circular 
cone  will  be  used  in  this  work.  It  is  suitable  for  describing  influence  functions  since  it 
tapers  off  at  the  radius  of  influence  and  has  a  peak  at  the  center  of  the  actuator.  An 
influence  function  is  displayed  in  Figure  87,  where  Ir  represents  the  radius  of  influence. 


INFLUENCE  FUNCTION 
z 


x 

Figure  87.  Influence  function. 


A  functional  description  of  the  influence  function  is  as  follows: 

Ik(X,Y,XA,YA)  =  [l--^V(X-XA)2+(Y-YA)2l  u[lf-{(X-XA)2+(Y-YA)2}]  (246) 
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where  X,  Y  =  spatial  location  on  the  mirror,  XA,  YA  =  spatial  location  of  an  actuator,  and 
u(«)  is  a  step  function  which  specifies  the  domain  of  each  influence  function. 

In  this  work,  Ir  has  been  chosen  equal  to  the  actuator  spacing,  which  implies  that 
the  actuators  are  independent,  i.e.,  the  actuators  only  effect  their  particular  mirror  segment 
and  do  not  effect  neighboring  segments.  An  actuator  map  with  numbered  actuators  is 
shown  in  Figure  88.  Actuators  1,  7,  43,  and  49  are  inert  while  all  other  actuators  are 
active. 


Y  (cm)  ACTUATDR  MAP 


Figure  88.  Actuator  map. 
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Considering  each  term  in  the  vector  m2T,  the  following  integral  is  generated  by 


substituting  eq.  (246)  into  eq.  (245) 


m, 


dxdy 


(247) 


where  ni2j  denotes  the  jth  term  in  the  vector  iri2T  and  ly,  uy,  lx,  ux  are  the  upper  and  lower 
limits  for  the  y  and  x  integrations  and  are  determined  by  the  step  function  in  eq.  (246). 
Performing  the  integration  in  eq.  (247)  yields 


[4-il]K-iy]+^f-Ki.+3xAi;+i;]. 


m2i  =- 

2j  tuR1 


3RI 


(248) 


In  the  determination  of  eq.  (248),  the  integrals  in  eq.  (247)  were  converted  to  equivalent 
polar  integrals.  In  doing  so,  the  polar  radius  was  allowed  to  vary  with  the  spatial  location 
X.  Utilizing  this  fact  along  with  the  actuator  map  in  Figure  88,  it  can  be  determined  that 
all  rows  of  actuators  will  have  the  same  nrbj  terms.  Therefore,  using  an  actuator  spacing 
of  0. 1  cm  gives 

mLPi  =  [0.08909  0.02339  0.00169  0  -0.00569  -0.03933  -0.12509] 


m 


2  rep  2 


2=[0  0.02339  0.00169  0  -0.00569  -0.03933  0] 


(249) 


m2  _  J^[m2rep2  m2repl  m2repl  m2repl  m2rep1  m2repl  m2rep2]  (250) 

where  m]repl  and  m]rep2  represent  the  repeated  entries  in  the  m2T  row  vector.  Using  the 
same  type  of  analysis,  m3T  becomes 
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0 

0.02339 

0.00169 

0 

-0.00569 

-0.03933 

0 


0.08909 

0.02339 

0.00169 

0 

-0.00569 

-0.03933 

-0.12509 


(251) 


B.3  Output  Equation  Matrix  N 

tT 

Recall  from  the  work  in  Chapter  5  that  N  =  Pxl  Px2  •  Pxp  Py,  PY2  •  Pyp  • 

Therefore,  the  vectors  pxj ,  pyj ,  j  =  1,  2,  ...,  p,  must  be  determined,  each  of  which 

contain  two  entries,  one  for  x-tilt  and  one  for  y-tilt.  In  the  case  of  x-  and  y-  tilts  only,  the 
N  matrix  becomes 

Pi,  (!)  Pxi  (2) 

PL  (0  PL(2) 

Pip  (1)  Pip  (2) 

Py,  CD  Pyl  (2) 

Py2  (1)  Py2  (2) 

Pyp  CD  Pyp  (2) 

where  pxkT(  1 )  and  pxkT(2)  are  the  first  and  second  elements  in  the  pxkT  row  vector.  From 
eq.  (157), 

9 


Wx  (X  -  Xs ,  Y  -  Ys )  -  Z^(X)  dXdY ,  n  =  1,2.  (254) 

x‘”~ 
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Taking  the  partial  derivatives  of  the  Zemike  polynomials  yields 


2x  az2  2  az2 

2  =  —  =>  -r-^-  =  —  ,^—L  =  0 

2  r  ax  r  aY 
_  2Y  az3  _  az3  _  2 
3“  r  ^  ax  _u’ aY  —  r 


(255) 


Therefore,  it  is  easily  seen  that  the  first  entry  in  each  p  ^  row  vector  is  zero,  while  the 
second  entry  in  each  i  is  zero.  In  terms  of  the  other  entries,  consider  first  the  row 
vectors  pY  .  The  integral  in  this  case  becomes 


ff  1 

2  In  2-1 

f,  2(Y-Ys)^ 

In 

f.  2(Y-Y,)'| 

■JJ  2 In 2 

Va  J 

l  Va  ; 

f 

1  + 
V 


2iY-Y^ 

Va  j 


In 


1  + 


2(Y-Ys) 

Va 


—  dXdY. 
R 


(256) 


Performing  the  integrations  in  eq.  (256)  produces 


Pi(2)  =  - 


2A 

R 


1- 


(ln  2-1) 
In  2 


2  In  2 


(257) 


It  can  also  be  seen  that  (1)  will  be  exactly  the  same  as  the  expression  in  eq.  (257). 


Hence,  the  N  matrix  becomes 
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N  = 


2A 


R 

2A 


R 

2A 


R 


0 

0 

0 


pxl 


(In  2-1) 

1 

In  2 

2  In  2 

(In  2-1) 

1 

In  2 

2  In  2 

(In  2-1) 

1 

In  2 

2  In  2 

2A 

R 

2A 

R 

2A 

~R 


1- 


1- 


1 


(In  2-1) 

In  2  2  In  2 

(in  2-  l)  1 


In  2 


2  In  2 
1 


In  2  2  In  2 


pxl 


0 

0 

0 


pxl 


pxl 
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B.4  Derivation  of  Corrected  Wavefront  RMS  Phase  Distortion 

In  this  section,  it  is  desired  to  derive  the  expression  displayed  in  eq.  (212).  Recall 
that  the  RMS  phase  distortion  can  be  expressed  as 

^RMs(0  =  Vw(x)O2(x,t)dx.  (259) 

For  the  corrected  wavefront,  the  expression  for  the  phase  becomes 

^(corrected)  (-,t)  =  £_a<«»)(t)Zj  (X)  +  ^ &("0  (t)Zk  (x)  (260) 

j=2  k=2 

where  the  negative  sign  arises  because  the  conjugate  phase  is  being  applied  to  the  mirror. 
Note  that  both  summations  start  at  index  two.  This  is  because  Zemike  number  one  is 
piston  and  is  not  included  in  any  calculations.  Also,  only  modes  two  and  three,  i.e.,  x- 
and  y-tilts,  of  the  mirror  coefficients  are  included.  Even  though  an  accurate 
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representation  of  the  outbound  wavefront  requires  n  Zemike  modes,  only  tilts  are 
corrected.  Substituting  eq.  (260)  into  eq.  (259)  gives 


*'"’(•)  =  J|w(x) 


n+1 


-|2 


j=2 


k=2 


dx 


(261) 


Expanding  the  summations  produces 


0"d)(t)  =  -J  j  W(x)[-a<">  ( t)Z,  (x)  -  a‘0)  (t)Z3  (x)-  •  -a^,  (t)ZB+1  (x)  +  a<"°  (t)Z2  (x)  +  a<m>  (t)Z3  (x)fdx 


(262) 


Expanding  the  term  in  brackets  in  eq.  (262)  yields 


[•]  =  {a^0)}2Z2Z2  +af  a‘0)Z2Z3  +---+a'°,a(n°+)1Z2Zn+,  - a‘0)a‘m)Z2Z2  -a(20)a'm)Z2Z3 
+a<°>a<o)Z3Z2  +{a'0)}2Z3Z3  +-+a(3°)a(n0+)IZ3Zn+1  -  a(30)a<m,Z3Z2  - a(30)a<m)Z3Z3 

+  -+al0+>a<0,Zn+1Z2  +a(n°+)1a(30)Zn+1Z3  +•••+  {a(n°+>}2Zn+1Zn+1  - a^a<m)Zn+IZ2  -a<Mm)Zn+tZ3 

— 7  7 _ (m)  (o)  7  7  -! -In(m)\27  7  A -n(m)n(m)7  7 

a2  a2  Z,2Z,2  a2  a3  >^2^3  a2  an+lZ'2Zjn+l  +\a2  J  Z"2Zj2+a2  d3  ^2^3 

-a3m)a20)Z3Z2  - a3m)a30)Z3Z3 - afV^Z.Z^,  +  a<m)a<m)Z3Z2  +  {a<m)}2Z3Z3  (263) 

where  Zj,  j  =  2,  3, ...,  n+1  is  Zj(x) .  The  vector  dependence  has  been  dropped  to  simplify 

notation.  Now,  the  orthonormality  property  of  the  Zernike  polynomials  can  be  used. 
Recall  that 


|w(x)Zj(x)Zk(x)dx 


[1  if  j  =  k 
[0  if  j*k. 


(264) 


Substituting  eq.  (263)  into  eq.  (262)  and  performing  the  integration  using  the  property  in 
eq. (264)  gives 
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*“,(l)  =  A/h;",(,)  +  ar(.)f  +[-a'«(t)  +  ar)(t)f  +{ar(t)}!  +-  +  {aS,(0}’  (2«) 

which  is  the  result  shown  in  eq.  (212)  of  Chapter  7. 
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employment  of  advanced  control  and  estimation  concepts.  The  control  objective  is  to  apply  the  estimated  conjugate  phase  to  the  deformable  mirror  so  that,  at  the 
target,  the  outbound  wavefront  distortion  is  minimized  and  the  Strehl  ratio  is  maximized. 
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